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FOREWORD 

This document is Volume I of a two --volume report describing the i 
Reacting Multi-Phase (RAMP) Computer Code developed by the Advanced^ / 

- ' ' ' % ■ 5 '/' 

Technology Section of Lockheed's Huntsville Research & Engineering i, 

' ' ' ■ \ 'V''" 

Center. Volume I addresses the theory and numerical solution roS;.th€5'^ 

Gomputer code. Volume II describes the computer code along with the 
program input and output. 

Documentation of the computer code was prepared in partial fulfill- 
ment of contract requirements (Contract NAS9-14517) with the NASA- 
Johnsonv Space Flight Center, Houston, Texas, in support of Space Shuttle 
related exhaust plume applications . The contracting officer's technical 
representative for this study was Mr. Barney B. Roberts of the Aero- 
dynamics System Analysis Section. 

The authors acknowledge the efforts of a number of individuals, *^ho 
/ - contributed to the development of the RAMP code. These include Dr, Terry 
F. Greenwood and Mr. David C. Seymour of the NASA-Marshall Space Flight 
Center; and Messrs . Robert J. Prozan, Jon A. Freeman, L.R. Baker and 
A. W. Ratliff of Lockheed-Huntsville. Ideas and suggestions for improve- 
ment of the analysis are reflected by frequent consultation with these 
individuals. 

Companion documents to this report include a user ‘s manual for the 
RAMP code; a report which describes the modifications made to the NASA- 
Lewis TRAN 72 computer code; and documentation of a one -dimens ional 
solution which provides a supersonic startline for the RAMP code. These 
reports are, respectively: 
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"Supersonic Flow of Chemically Reacting Gas- Particle 
Mixtures"— Volume H RAMP — A Computer Code for 
Analysis of Chemically Reacting Gas -Particle Flows," 
LMSC-HREC TR P496555-II. 

"User's Guide for TRAN72 Computer Code Modified for 
Use with RAMP and VOFMOC Flowfield Codes," LMSC- 
HREC TM D390409.^^^ ^ 

"General One -Dimensional Flow of a Gas-Particle 
System," LMSC-HREC TM D390876. 
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SUMMARY 

' • ■ ■ . V-''” ' ■ ' ■ 

This report describes a numerical solution for chemically reacting 
supersonic gas -particle flows in rocket nozzles and exhaust plumes. The 
gas -particle flow solution is fully coupled in that the effects of particle drag 
and heat transfer between the gas and particle phases is treated. Gas and 
particulates exchange momentum via the drag exerted on the gas by the par- 
tides. Energy is exchanged between the phases via heat-''transfer (convec- 
tion and/or radiation between the gas -particle phases). 

Basic assumptions made in the development of the governing equations 
are similar to those employed by previous investigators. The primary ex- 
ception is the treatment of chemical effects in the gas phase. Thermo- 
chemistry calculations (chemical equilibrium, frozen or chemical kinetics) 
are shown to be uncoupled from the flow solution and, as such, can be soivjed 
separately. 

The solution to the set of governing equations is obtained by utilizing 
the method of characteristics . The equa^i ons cast in characteristic form 
are shown to be for ideal, frozen, chemical equilibrium 

and chemical non -equilibrium reacting gas mixtures. The characteristic 
directions for the gas -particle system are found to be the conventional gas 
Mach lines, the gas streamlines and the particle streamlines. 

The basic mesh construction for the flow solution is along streamlines 
and normals to the streamlines for axisymmetric or two-dimensional flow. 
The analysis gives detailed information of the supersonic flow and provides 
for a continuous solution of the nozzle and exhaust plume flow fields. Bound- 
ary conditions for the flow solution are either the nozzle wall or the exhaust 
plume boundary. 
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The particle distribution is represented in the numerical solution by 
a finite distribution of particle sizes. The particle limiting streamline con 
cept is utilized to define the region of influence for a particular particle. 
Particle physical and thermodj^namic properties are defined by the particle 
mass density, distribution of particle sizes and thermodynamic data. 

Presented is the development of the set of governing partial differ- 
ential equations for the gas -particle^ system. The governing equations are 
cast in characteristic form and the corresponding difference equations 
formulated. The numerical solutions for the various point types are de- 
scribed and the corresponding steps of each solution outlined. 


Vi : 
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Symbol 

A 


a 



dv 


A 

e 

E 

F ■ 

f 

G 

H 

k 

1 

sp 

K 

k 


SYMBOLS AND NOTATION 
Description 

parameter defined by Eq. (2.46), 1/sec; area, ft 
speed of sound, ft/sec 
parameter defined by Eq. (2.86c) 

■ii ^ ^ 

paranrieter defined by Eq. (2.83), ft /sec /°R 

drag coefficient, dimensionless 

specific heat at constant pressure, ft^/sec^/°R 

specific heat at constant volume, ft^/sec^/°R 

drag force, Ibf 

substantial derivative 

an element of volume, ft 
unit vector, dimensionless 

energy in the gas -particle system control volume, 
ft^/sec2 

interpolation factor, dimensionless; or total force 
acting on the system defined by Eq. (2,30), Ibf 

drag coefficient parameter ) V r 

dimensipnless * Stokes 

Nucselt number parameter defined by Eq. (2,77), 
dimensionless 

2 2 

total enthalpy, ft /sec 
local enthalpy, ft' /sec^ 
specific impulse, Ibf-sec/lbm 

parti cle keat-=t^an sfer film coef fici ent defined by 
Eq. (2.76), lbm/seG/9R/ft 

thermal conductivity of gas, lbm/sec/*^R 



PREGBDING PAGE BLAHE K05' PI^^ 
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Symbol 



NG 

NP 


NS \ 

Nu 
n 

,'P' 

p. 

I 

Pr 
Q 

q 

R 
Re 

rj 

S 

T , 
t 

U, V 

V 
W 



Description 

Mach number, dimensionless; or momentum, lbm-£t/sec 
mass flow rate, slug/sec 


density of a partial 
particles, slug/ft 


e based 


on a unit volume of 


index denoting number of discrete gaseous species 
considered 


index denoting number of discrete particles considered 


index denoting number of,discrete species 
considered = NG+NP 

Nusselt number, dimensionless 

unit normal vector 

2 . 

pressure, Ibf/ft 

2 

partial pressure of species i, Ibf/ft 
Prandtl number, dimensionless 

heat transferred to or from a gas -particle sy stem 
control volume, Btu 

velocity, ft/sec 

gas constant, ft^/sec^/°R 

Reynolds number, dimensionless 

particle radius , ft 

surface area of a gas -particle system control 
volume; ft 2; entropy, ft^/sec^/op 

temperature, °R 

time, sec 

velocity components, ft/sec 

3 

volume of a gas -particle system control volume, ft 

work performed on or by the gas-particle system 
control Volume, ft^/sec^ 

xLL 
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Symbol 


Description 


w 

X, r 
Greek 
a 

P 

r 

<p 

6 

r 

V 

0 

u 

X 

P 

a 

cr 

L 

X 

8 


production, rate, Ibm/sec 
position coordinates, ft 


accommodation coefficient, dimensionless; or Mach 
angle;^ radians 

characteristic slope, radians 
specific heat ratio, dimensionless 
any extensive quantity of the element dv 
Kronecker delta 

0, 1 for two-dimensional or axisymmetric flow, respectively 


2 

viscous stress tensor, Ibf/ft 
emissivity’, dimensionless 
nabla 

inclination of the flow vector \ 5 ^ith respect to the 
x-axis, radians 

chemical potential, cai/gm 

2 

viscosity, Ibf-sec/ft 

equation modifier 
2 

density, slug/ft 

2 

surface stress tensor, Ibf/ft 

2 3 

Stefan -Boltzman constant, ft /sec 
particle stream function, Ibm-sec/ft 


indicates summation. 


N 


j = l 


species mass fraction 
indicates partial derivative 


xiii 
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Subscripts 

I 

i 

L 

m 

n 

w 

X, V 

Superscripts 


Description 

indicates initial data surface 

indicates the quantity pertaining to species i 

denotes local surface conditions 

indicates the quantity pertaining to the gas -particle 
mixture 

indicates n^^^ data surface 
denotes nozzle wall conditions 

denotes partial differentiation in the x and y (or r) 
directions 

denotes a vector quantity 

denotes an average value over a step length 
indicates the quantity pertaining to a particle specie 


xlv 
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Section 1 
INTRODUCTION 


Solid propellant motors are frequently utilized to provide launch vehicle 
boost and stage separation. Most solid rocket motor propellants contain metal 
additives which increase the energy content of the s and suppress com- 

bustion pressure instabilities. The presence of metal additives results in 
condensed products in the nozzle -exhaust flow fields. This can create a 
number of adverse effects. 

Condensed products (i.e., particulate aluminum oxide, etc.) being inert 
can do no expansion work. Consequently the particles are accelerated through, 
the flow field via drag exerted on the gas by the viscous shearing X^tiono be- 
tween the gas and particulate phases. Since the particles do no 
work, the gas phase cools more rapidly than the particles. The particles at 
any given location are thus at a higher temperature than the gas so that heat 
is thus transferred from the particles to the gas by conduction and radiation. 
The net result is that the gas phase expends useful work in accelerating the 
particles while acquiring heat from the particles. This is an irreversible 
non-adiabatic process in the gas phase expansion. 

Exhaust plumes create hostile environments to surfaces immersed in 
the flow. These are in the form of structural loads, contamination and heat- 
ing. Heating results from both gas and condensates impinging on the surfaces. 
Particulates in the flow can also cause surface erosion and other structural 
damag e. 

Exhaust plume applications and plume impingement problems typically 
occur during launch, staging and rendezvous. The current Space Shuttle de- 
sign offers several illustrations of areas where exhaust plume flowfield prop- 
erties and flow structure must be known. This is schematically illustrated 

1-1 
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in Fig. 1-1 which indicates potential plume related problems from solid rocket 
motors. During the Space Shuttle launch the solid rocket booster motor ex- 
haust plumes impinge on the launch stand hardware. While the launch vehicle 
is in the vicinity of the launch pad, reradiation from the launch pad to the 
orbiter is a potential problem. The solid rocket motor exhaust plume affects 
the shuttle base drag and aerodynamics during the boost phase. The solid 
rocket separation motors mounted on the booster motor’s forward and aft 
ends used to effect stage separation can potentially subject the orbiter and 
external hydrogen/oxygen tank to a rather severe environment. 

/(( 

Exhaust impingement applications require a detailed definition of the 
gaseous plume structure as well as particle trajectories and dynamic prop- 
erties. Plume simulation studies require that plume shapes be known. The 
exit plane pressure must be known along with the gas thermodynamics. To 
adequately describe the flowfield properties, the mutual effect of gas on 
particle and particle on gas must be calculated. Also ’’real gas” effects can 
be significant and should be included in the gasdynamic considerations . 

Solutions for the supersonic flow of a gas -particle mixture follow one ^ 
of two approaches . These are: (a) a fully coupled solution in which momentum 
and energy are exchanged between the gas and particle phases (Refs, 1, 2 and 3), 
and (b) an uncoupled solution (Ref. 4) in which particle trajectories are traced 
through nozzle- exhaust plume flows. The uncoupled solution considers ’’real 
gas” effects; however, it treats only gas effects on the particle and is more 
applicable for low aluminum content propellants. The coupled solutions are 
primarily performance oriented. These solutions readily permit treatment 
of highly aluminized propellants but are restricted to constant thermodynamic 
properties. The performance code developed by Kliegel (Ref. 3) was subse- 
quently utilized to trace liquid droplets in predicting contamination to surfaces 
(Ref. 5). 

Applications discussed previously require a knowledge of the flowfield 
structure. Previous studies (Refs. 6, 7 and 8) indicated the need to include 
the treatment of ’’real gas” effects in nozzle -plume calculations. The 

'■ 1 - 2 , 
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decision was subsequently made to extend the coupled gas-particle solution 
(Refs. 1 and 2) to include the treatment of gas thermochemistry (chemical 
equilibrium and chemical kinetics) . The choice of computation scheme and 
computer code becomes more obvious after examination of anticipated 
applications. 


The method-of-characteristics has proven to be a reliable numerical 
solution for nozzle -plume applications (Ref. 9). Existing gas -particle per- 
formance codes utilize the method of characteristics which indicated that an 
operational nozzle-exhaust plume code could be obtained with a minimum of 
numerical solution development. The gas -particle capability has been in- 
corporated into a streamline -normal code (Ref, 10) which utilizes the method 
of characteristics to solve for local flow properties . There are several 
reasons for choosing this approach. These include: 

• The multiple shockwave capability is treated quite 
readily. 

• Numerical difficulties encountered in highly expanded 
flow are circumvented. 

anism for tracing particle streamlines. 

• Transition flow between the continuum and free molecular 
flow regime is more readily treated. 

The code in its present form will handle the flow problems which exhibit 
or have any of the following characteristics: 


• Supersonic In viscid 

• Highly Underexpanded\ioz 2 s\Vs 

• Highly Overexpanded Nozzles 

• Shock Waves * 

• Sliplines. 

• Solid ^ 

• Pressure Boundary 

• Nonequilibrium Chemistry (Finite Rate) 
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• Equilibrium Chemistry 

• Ideal Gherni|fry " 

• Free Molecular Flow 


• Two Phase Flow 

• Oxidizer-^ to-Fuel (o/F) Ratio 
Gradients 

' 


This report presents a detailed development of the equations governing 
the supersonic flow of a chemically reacting gas -particle mixture. Develop- 
ment of the governing relations follows to a large extent the work of Kliegel. 
Basic assumptions are stated and the set of partial differential equations de- 
scribing the gas -particle system are developed. These relations are then 
cast in characteristic*' form and the corresponding difference equations 
written. Details of the nurnerical solution are described for the various data 
point types. The presentation is then concluded with a description of the 
numerical integration of the conservation equations. 


In addition to a description of the above analysis techniques/ appendixes 
are included which discuss; (1 ) particle drag and heat transfer coefficients ; 
(2) non -isoenergetie gas-phase flow treatment; (3) chemical equilibrium calcu 
lations in gas-particle flows; (4) non-continuum flow expansions; and (5) inte- 
gration of the finite -rate chemical kinetic equations. 
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Section 2 

FUNDAMENTAL EQUATIONS FOR STEADY FLOW 
OF REACTING GAS-PARTICLE MIXTURES 

2.1 BASIC ASSUMPTIONS FOR THE GOVERNING EQUATIONS 

The flow of a gas -particle mixture is described by the equations for 
conservation of mass, conservation of momentum and conservation of energy. 
In the gaseous phase the state variables, P, p, R and T are related by the 
equation of state while for the particulate phase the equations are for the 
particle drag, particle heat balance and the particle equation of state. De- 
velopment of these equations is based on the following assumptions: 

1. The particles are spherical in shape. 

2. The particle internal temperature is uniform. 

3. The gas and particles exchange thermal energy by con- 
vection and radiation (optional), 

4. The gas obeys the perfect gas law and is either frozen 
and/or in chemical equilibrium, or is in chemical non - 
equilibrium. 

5. The forces acting on the control volume are the pressure 
of the gas and the drag of the particles . 

6. The gas is inviscid except for the drag it exerts on the 
particles . 

7. There are no particle interactions. 

8. The volume occupied by the particles is negligible. 

9. There is no mass exchange between the phases. 

10. A discrete number of parti^ each of different size or 
chemical species, is chosen to represent the actual con- 
tinuous particle di stribution. 

11. The particles are inert. 

These are basically the same as originally stated by Kliegel except for the 
provision to: (a) calculate the radiation exchangd between the gas and particle 
phase (assumption 3); (b) the gas phase can be in either frozen,; in chemical 
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equilibrium, or chemical non -equilibrium (assumption 4). The toal mass 
and energy of the gas-particle system are not constant since provisions are 
made for particle streamlines to penetrate the exhaust plume boundary. 

2.2 GOVERNING EQUATIONS FOR THE GAS-PARTICLE MIXTURES 
2.2.1 Continuity Equation 

- If ' 

The various forms of the continuity equation for chemically reacting 
gas -particle mixtures are derived in this section. The forms of the equation 
derived, in order are: 


1. Species Continuity Equation (2.6) 

2. Global Continuity Equation /(Steady State) (2. 13)/(2.27) 

3. Gas Continuity Equation' (2.19) 

4. Particle Continuity Equation (2.20) 

5. Species Continuity Equations as a Function of (2.25)/(2.29) 

Mass Fraction/(Steady State) 

6. Particle Continuity Equation for Each Particle (2.26)/(2.28) 

species j/(Steady state) 


Consider a chemically reacting gas -particle mixture flowing through 
some arbitrary stationary control volume, v, (Fig. 2-1) bounded by the con- 
trol surface, S. 



Fig. 2-1> Control Volume for the System of Equations 
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For such a flow system, Reynold's Transport theorem states that; 

^ / AdV = / ^dVt /A^f.SdS; (2.1) 


V 


where A = any arbitrary parameter. 


Applying the Gauss, or Divergence theorem. 

J" 'S> • n dS = ^ V . B dV ; 
S V 

where B = any arbitrary vector quantity 


to Eq. (2. 1) yields; 


& 


m / A dV = d'Y + f V . A^ dV . 


( 2 . 2 ) 


(2.3) 


Furthermore, the substantial or Euler 's derivative may be written in the 
form; 

DA 3A . 

Dt 8t * 


(2.4) 


In a flow system of gas -particle mixtures in which chemical reactions 
take place, the principle of conservation of mass of each chemical species 
may be written as; 


^ y*p^ dV - y* w. dV = 0 i = l.NS; 
V V 


(2.5) 


where - production rate of species i du^ to either Internal or external 
sources such as chemical reactions and mass additions. / 


Applying Eq. (2.3) to the first term of Eq. (2^$) 


^ A 


dV = 


Ei 

8t 


dV + 


/ ^ • Pi^i 


dV i = 1,NS , 


V 
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substituting the result back into Eq. (2.5) 

9P; 


J' dV +y*- V • ^ w^dV = 0 i = 1,NS , 


V 


and rearranging terms yields 


/ ^ = - / 
V 

" TT 


dV = - I V • p. q dV + 

V V 




dV 


i = 1,NS 


■ The above equation simply states that, for species i, the mass rate of 
increase inside the control volume is equal to the net rate of mass flow into 
the control volume plus the rate of mass produced due to chemical reactions 
and mass addition. In the following, we shall, however, exclude mass addition 
from external sources. 

Combining the integrands under one integral 

y* +V.P.5; -wAdV = 0 i = l,NS, 

V ' 

and noting that the volume V iinder consideration is arbitrary; the only way 
for the above equation to be valid for all V is for the integrand to vanish. We 
therefore have 


9pi . 

-^ + V Pi q. - Wi 


= 0 i = 1,NS . 


(2.6) 


Equation (2.6) is known as the species continuity equation. It is valid 
for each chemical species at each internal quantum state. We shall, however, 
assume that the various internal modes of motion are fully excited and are 
in equilibrium with each other . It is well known that this is approximately 
the case for the translational and rotational degrees of freedom where the 
equilibrium value is attained in a few collisions. In general, the vibrational 
degr je of freedom approaches the equilibrium state somewhat more slowly, 


2-4 


LOCKHEED - HUNTSVILLE RESEARCH 4 ENGINEERING CENTER 


except at very high temperatures. As chemical reactions usually occur at 
high temperatures, this approximation is often justified. The global continuity 
equation can now be derived by summing Eq. (2.6) for all species present. 


or 



( 2 . 7 ) 


To arrive at the final form of the global continuity equation, it is necessary ' 
to take a closer look at the flow system. 


In the flow system analyzed, the time variation of the thermodynamic 
functions is slow compared to the longest relaxation time of the system. 
Therefore, the assumption that thermodynamic local equilibrium, exists can 
be made. In such a system, the thermodynamic quantities for the non equi- 
librium system are the same functions of the lodal state variables as the 
corresponding equilibrium thermodynamic quantities. Therefore, the partial 
specific quantity in an arbitrary volume element dV of a nonequilibrium 
system may be defined by the equilibrium relation 



The quantities being held constant during the differentiation are the locally 
defined temperature T, pressure P, and the masses m^ of the other NS-l 
species . 


The specific quantity (4> per unit mass) is given in terms of the 


partial specific quantities by 


NS 


p ^ = y^p. f i = 1,NS; 

^m ’^m / .. 1 I 
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REPRODUCIBILITY OE- THE 
ORIGINAL PAGE IS _ 

(6/^ 


where the relationship between the mass densit/ of the mixture and the 


species density is given by 


m 


NS 


. = fluid density 


( 2 . 8 ) 


i=l 


The specific velocity, q , specific enthalpy, h , and specific internal 

m ^ 


energy, are given by; 


1? 


and 


NS 

^m ^m ~ *^i 

i=l 

NS 

p h = p. h. 

^m m 1 

i=l 


m 


P e = p. e. 

*^m m iJL# 1 


NS 

E 

i=l 




(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


Substituting Eqs. (2.8) and (2.9) into (2.7), 


9P. 


NS 


m 


at 


V • p q - w. = 0 , 
*^m ^m Z-/ 1 

i=l 


and noting that for a closed system, the total mass for chemical reactions 
is conserved 

NS 


i=l 


w. = 0 ; 
1 


( 2 . 12 )' 


Equation (2.7) reduces to 


9P 

■# *^m V = ° 


(2.13) 


Equation (2.13) is known as the global continuity equation. It can be written 
alternatively as • , 

^ V = 0 


by means of the substantial derivative. 
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Equations (2.8), (2.9) stnd (2.12) can be rewritten in terms of gaseous 
and particle species as follows: 

Equation (2.8) 


NS 

"m = Z Pl 

i=l 

Equation (2.9): 

NS 

i=l 

Equation (2.12): 

NS 

“ 

i=l 


NP . 

(2.14) 

j=l 


w . . 

= pq+2jp^q^ (2.15) 

j=i 


NP' 

+ 2 = 0 (2.16) 
j=l 


Assuming that reactions of the form A(gas) + B{gas) ‘TZSi C(particle) + D(particle) 
do not contribute significantly to the system (interchange of phases is negligible 
due to chemical reactions), Eq. (2.16) implies that 

NP 

w = 0 and = 0 . (2.17) 

j=l 

Substituting Eqs. (2.14) through (2.17) into Eq. (2.7) 

0 0 



and rearranging terms yields 


NP 


|& + V Pq +2^^ +V.pj?) = 0 


(2.18) 


2-7 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 








Since the gaseous species and particle species do not interchange phases, 
both contributions must vanish. Therefore 


+ V • p q = 0 (2.19) 

and 

=0 . ( 2 . 20 ) 

j=l 

Equation (2.19) is known as the gas continuity equation and is valid when 
applied to the system of gas species. Equation (2.20) is known as the particle 
continuity equation and is valid when applied to the system of particle species. 


The species continuity Eq. (2.6) may be written in a more convenient 
form by making use of the gas continuity equation and introducing the mass 
fraction of species i. In general, the mass fraction of species i may be 
defined as; 



( 2 . 21 ) 


Substituting X.^ into the species continuity Eq. (2.6) 




m 


) +V 


Pm 4i 


W." 

1 


= 0 ; 


1,NS 


expanding the first term 

9p 

1 * m 

P “aT” “51“ +V * X. p q. - w. 
^m 9t '1 9t 1 '^m n i 


= 0 i = 1,NS 


(2.21a) 


and applying the vector identity 


^ • ^i ^m ‘li = ^ • (^i<^m ^i» = ^^i • ^m ^i + ^i * ^m ^i^ 
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i = 1,NS 


( 2 . 22 ) 


results in 


rax, 

i 9t 




+ V • p„ q 
m 


i) - »i ' » 


\v. i; 


By limiting the system to only gas phase reactions (no reactions in the particulat* 
phase) the species mass fraction may be redefined as: 


= p/p 

Assuming all gaseous species have velocity q^, = q, Eq. (2.22) becomes 


(2.23) 




i = l,NS (2.24) 


Applying the gas continuity equation to the above expression the second term 
is set equal to zero and Eq. (2.24) becomes: 




i = l.NS 


The final form of Eq. (2.24) is obtained by applying substantial derivative, 
Eq. (2.4) to yield 


DX. 

p - = 0 i = I, NS ■ species continuity 

equation 


(2.25) 


Under the above assumption that there are no particl'6 reactions, Eq. (2.20) 
for each particle species j becomes; 


4 V •p'^q^ = 0 j = 1,NP particle continuity 

equation 


(2.26) 
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iT 


For steady state flow processes, Eqs.(2.13), (2.26) and (2.25) become: 
Equation (2.13) ^ 


= 0 


global continuity- 
equation 


(2.27) 


Equation (2.26) 


V • = 0 j = 1,NP particle continuity 


(2.28) 


equation, and 


0 i = 1,NS species continuity 
equation 


(2.29) 


Equation (2.25) 

p J. . VA'. - w. 

O ' 

2.2.2 Momentum Equation 


The various forms of the global momentum equation for a chemically 
reacting;gas -particle mixture are derived in this section. The forms of the 
equation derived, in order are: 


1. 

General Global Momentum Equation 

(2.42) 

2. 

General Global Momentum Equation 
with particle drag effects 

(2.49) 

3. 

Global Momentum Equation for steady 
state, inviscid, no body force flow 

(2.50) 


The linear momentum of an element of mass m is a vector quantity 
defined as mq . The fundamental statement of Newton's law for an inertial 
reference is given in terms of momentum as 

(2.30) 


D 






r:T 




’ )• 


I 




ill 


Si 




wr 




F - (m q ) = 0 
Dt ' m ^m 


To derive the general global momentum equation the two terms of Eq. 
(2.30) will be analyzed separately. 
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The total force acting on the system may be represented as the sum 
of the surface force distributions (force distributions acting on the boundary 
of the system) and body force distributions (force distributions acting on the 
material inside the system) as follows; 


= j <j •nds +y p. f. 

V i=l 


dV 


(2,31) 


The surface force term may be converted into a volume integral by applying 


Gauss' theorem, Eq.(2.2). Therefore, 

r = /. NS 

F = J V . cr dV + / p. f. dV 

V i=l 


(2.32) 


The surface stress tensor, a may be written in the form: 


and 


a = -P6 + T 


V • <7 = - VP + V . T 


(2.33) 


(2J34) 


Substituting Eq. (2.34) into Eq, (2.32), the total force acting on the system 
may be written in the form 

NS 

F = y* (- VP + V . T) dV + y 23 ^i h ’ 


or, upon combining terms 

F 


-f (- 

V ' 


VP + V , r + 


V i=l 


NS ' \ 
i = l / 


dV 


(2.35) 
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The momentum term m q in Eq. (2.30) may be rewritten in the 

xn m 


form 


^ NS 

m q = r p. q. 

m^m J 


: q. dV ; 


V i=l 


from which the substantial derivative can then be written 


' Dt 

W 

'A 


D D /* 

Dt <“mV> = D? y £ "i % • 


NS 

£ 

V i = l 

Applying Reynold's transport theorem Eq. (2.3) to Eq. (2.36) 

NS _ NS 

£ 

i-1 

Y V 

\ 

rearranging terms 

NS NS 


(2.36) 


^ .ft NS NS 

•”mV> = / 8t£<>i’i ; 


JJt m m 


' / It £pi^i+^‘£'>i'^«i( 

' i = l i=l 


(2.37) 


and rewriting the results in terms of the gaseous and particle species results 

I 

in: ■ 


ft 


= f jif (p q + V • (p 4 q + p^q^^jjdV 

V ' ' j=l C ' ' j=l ' 


or, 


Dt 


( il NPp ,v 

= / lt<p5+^-'AA^+£^'P^5^)+^-P'?^ <>V <2-38) 

V ' , : •/: j=r ' 
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Substituting Eqs, (2.35) and (2,38) into Eq. (2.30), dropping the integral notation 
and rearranging the terms yields 

NPr ft . . . . ..i _ NS 

^ (P^ ?) + V • pJ ? qJ = -VP +V . T + 2 (2.39) 

j=l^ i=l 


^(p q) • P q q +' 


Expanding the first two terms of Eq. (2.39) 


_8_ 

8t 


(pq) +v • Pqq = P^ +q|^ +pq* vq + qv *p q; (2.39a) 


and then applying the gas continuity Eq. (2.19) 


^ + V • P q = 0 


9t 


to the results, yields 


at 


(P q) + V • p q q = 


0 Q — i 

P (-^ + q • vq) 


(2.40) 


The third term of Eq. (2.39) may be manipulated in a like manner using the 
particle continuity Eq. (2.26) to yield; 

(P^ ^) + V • P^ q^ q^ j = + ? • Vq^) 

Substituting Eqs. (2.40) and (2.41) back into Eq. (2.39) results in; 

-’NP ^ _ NS 

p(^ +q •vt)+ 2!^pM-^ + ?j.v?n= -VP+V .T+ XPi'fi ; 

' ' j=l • ' i=l 


(2,41) 


or in terms of the substantial derivative, Eq. (2.4), 

p dT <1 + a E"! 'i 

j=i i=i 


(2.42) 
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Equation (2.42) is the general global momentum equation. In this equation, 
the five terms are, respectively; (1) the non -stationary and convective rate 
of change of momentum per unit volume of the gaseous species; (2) the ni^n- 
stationary and convective rate of change of momentum per unit volume of the 
particle species; (3) the net hydrostatic pressure force; (4) the, viscous stress 
force acting on the surface of the unit volume; and (5) the body forces per unit 
volume. 

Since it is assumed that all forces are negligible, except that of the 
pressure of the gas and the particle drag, the only force exerted on the 
particle is a drag force d| caused by the relative motion between the gas 
and the particle. Figure 2-2 illustrates a spherical particle in a gas flow field 



Fig , 2-2 - Particle Momentum and Heat Transfer Model 

To determine the effects of particle drag on the global momentum 
equation, Newton's Second Law is applied to the drag forces acting on 
particle species j. 


v : ; ■ 2-14 ; 
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or rewritten in terms of particle density 

In general, the drag force on a spherical particle may be written in the form: 

jr(r'^)^ (ip Aqj )Aq^|) 


where 


Aqj = q - . 


(2.43) 


Equating the two expressions for particle drag 

|-s-(rj) m'^j^qj = cj^ ^(r'^) ip Aqj jAq'^l; 


and rearranging terms yields: 

D -M 
Dtl 


3 /P C|)\ 

" « W rV 


a?J|a?| 


(2.44) 


Equation (2,44) can be simplified by referencing to the Stokes flow regime 
(Reynolds number is Iqssrthan 1) in the following manner. 


Define 




fj = 




Stokes 


where 




Stokes 


r k’ = 2 


Therefore, 


(-jj _ 2^fj^ - 24 fj 
- Re ~ ^ 


^2 r j p Ia^I ^ 


(2.45) 
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Substituting Eq. (2.45) into Eq. (2.44) and defining the parameter 


j _ 1 


A-* = 




. .2 


(2.46) 


the momentum balance becomes 


^qj = 


(2.47) 


Expanding the substantial derivative by means of Eq. (2.4), Eq. (2.47) becomes 


= aS A? = ^ . 


Dt 


(2.48) 


Substituting Eq. (2.47) into the general global momentum equation results in: 
^ NP _ NS 


J=1 ■ i|l 


(2.49) 


Equation (2.49) is the general global momentum equation drag 

effects. 

: / 

Fo:iL^/he be described as steady statef inviscid 

and no body forces present Eq. (2.49) becomes: 

NP. ^ . 

p q «Vq + E P'^ Aqj = - VP. (2.50) 

We note at this point that chemical reactions 'do not alter the forms of the " 
global continuity equation or equations of motion. 


Equation (2.50) will be expanded now for later use in the derivation of 
the energy equation. For an arbitrary vector A the following identity^ exists 


1 


A • VA = V (A I A) - Ax V x A) 


(2.51) 
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Applying this identity to the first term in Eq. (2.50) 

P q • Vq = 2 P V(q • q) - pq q) ; 

substituting the results back into Eq. (2.50) 

NP 

^P V(q »^) - P ^ x(v xq) + A? = - VP ; 

j=l 


and rearranging terms yields 

NP \ 

•| V(q . q) = q x(V xq) ~ p ^ ~ ^ 

j = l 


^ (2.52) 


2.2.3 Energy Equation 

The various forms of the energy equation for a chemically reacting 
gas-particle mixture are derived in this section. The forms of the equation 
derived are: 


1. General Global Energy Equation Neglecting the Effects of 
Radiation 

2. General Particle Energy Balance Equation 

3. General Global Energy Equation with Radiation Effects 



Global Energy Equation with Radiation Effects for Flow 
Described as Steady State,' Adiabatic, Inviscid, and No 
Body Forces Present 


(2.74) 

(2.84) 

(2>85) 

■r 

(2.105) 


The most general form of the energy equation, according to the first law of 
thermodynamics or the law of conservation of energy can be written as 


DE dQ dW 
Dt dt dt 


(2.53) 
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„The energy E in the control volume may be written in terms of the species 

'Mh -. - ' - ■ 

“-present as 

^NS' 


where 


V i=l 


qr q. 

e. = i. + —■ + Z. i = 1, NS 

1 1 2 X 


(2.54) 


i. = internal energy of species i per unit mass 

qf . V 

— = kinetic energy of species i per unit mass 


Z. = potential energy of species i per unit mass. 


8c ^ 


The substantial derivative can then be written 


DE 

Dt 


r» /• NS 

V i=l 


and expcxnded by applying Eq. (2.3) to yield 


DE 

Dt 


- NS NS 


V i=l 


i=l 


Recalling that the total force acting on the control volume according to Eq. 
(2.31) is 

,NS 


Pj ^ dv ; 


V i=l 


and applying the definition for the rate at which work is performed on or by 
the control volume dW/dt = -q • F, 
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one obtains 


NS 


^ = -y*? .<T .n ds .y* ^ . p. dv 

S V i=l 


The surface integral may be converted into a volume integral by applying 
Gauss ' theorem, Eq. (2.2). Therefore, 


dW 

dt 


NS 


= - y* V (a . q) dv - y ^ q. . p. t dV . 


i=l 


Recalling the definition of the surface stress tensor a , Eq. (2.33) 


cr = - P6 + T , 


and taking the dot product 


<j • q = - Pq + T • q , 


the first volume integral in the above equation may be written as 


yV •(- Pq + T «q) dV . 
V 


Substituting Eq. (2.56) back into the expression for dW/dt 

NS 


dW 

dt 


= - f V • (- Pq + T« q) - r ]^q. . p. t dV , 
V y i=i 


combining integrals 


N S ^ ^ 

^ = - /* I V • (•- Pq + T* 5 + 2^ Ai • Pi , 

-- i=l 


and expanding terms yields 


NS 


V 


dV 
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(2.56) 


(2.57) 






The rate at which heat is transferred to or from the control volume may be 
written in the form 

f =/3.ndS. 

S • ■ . ' 

' ex'" 

or in terms of a volume integral by applying Gauss' theorem, Eq.(2.2) 




^ = y* V «Q dV , (2.58) 

V 

where Q is the conduction heat flux vector. 

Substituting Eqs..(2.55), (2.57) and (2.58) into Eq. (2.53) and dropping the 
integral notation results in 

fl NS NS ^ ^ _ NS ^ 

~ + ^ *Pq - V (r .q) - \ = 0* (2.59) 

i=l i=l i=l 

Equation (2.59) is the general global energy equation neglecting the effects of 
radiation. This equation is based on a unit volume. The six terms are, re- 
spectively: (1) and (2) thetotal rate of increase of internal, kinetic and potential 
energy caused by local and convective changes; (3) the heat conducted to or 
from the control volume; (4) and (5) the work done on the fluid element due to 
surface forces; and (6) the work done on the element by the body forces. 

In deriving Eq.|,^?2.59), all that has been neglected is the energy tratts- 
ferred due to radiation, which if necessary, can be added to the equation as 
an extra term. Radiation effects will be determined later in this sectipn. 

The microscopic quantum effects (interchange of energy and mass)^are 
not considered in Eq. (2.59) as it has been implied by the classical first law 
of thermodynamic s. The general global, energy equation will now be expanded 
and written in terms of the gaseous and particle species present.^^^ ^^ ^^^ ^ 
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Recall that 


2 

q. q. 

e. = i. + — Z. i = 1,NS . 

1 I IL 1 

// 


If the potential energy term is; combined with the body force term, Eq. (2.54) 
may be rewritten in the form C 


e. = i. 

1 1 2 


i = 1,NS 


Defining the specific enthalpy of species i as 


h. = i. i = l.NS (2.61 

lip. 

where = partial pressure of species i 
and combining Eqs. (2.60) and (2.61) yields 

P. q? 

e. = h. - — + -Tj- 
1 i p. 2 

Next, define the total enthalpy per unit mass of species i as 

2 

H. = h. +-^ . (2.62 

1 1 2 

Therefore, the expression for becomes 

P. 

e. = H. i . 

X p. 

- . V. /r*' 

Substituting the expression for into Eq:'(2.59) results in 

NS j P.V / P.\ — 

-pT) -^]4i - V . Q tv . P^-V. (?-5 Pi'i = 


(2.62a) 
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fi 3 


/ 


The pressure P may be expressed as the sum of the partial pressures as 
follows 


NS 

p = • 

i=l 


(2.63) 


Substituting Eq. (2.63) into Eq. (2.62a) 


From the assumption that thermodynamic local equilibrium exists, the foilow- 
ing expressions may be written in terms of the gaseous and particle species: 

■ NS ^ . 

= PH+V;pJHJ (2.65) 


i=l 


j=l 


and, 


NS • 

2 ^ P ^ ^ 

i=l j=l 


( 2 . 66 ) 


JL 

at 


Substituting Eqs. (2.65) and (2.66) into Eq. (2.64) 
NP \ NP 


' '% 


iNJr ^ V IMP ^ ^ ^ ^ 

p H + --^ + V . ^pq H + hJ) - V. Q - V (r . q) 


\ 'll 


NS 

E^- pi^ " 0. 

i=l 
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II 




NS 

NS 

NS 


^Ep.h, 

-^ + ^-EpiqiHi- V- 

Pj ti - V . Q 

■‘Sr 

1=1 

i = l 

i=l 




NS 



+ V • Pq - V • (r • q) - 





i = l 


and noting that 

NS 


1 


v.EqiPi = v 

•Pq. 



i=l 



the above equation becomes 



^ NS 

NS 

_ ^ NS 

. Hi' 

lEPi «i - W + ^ -D - V -3 - V 

• (T • 5 - X) % * ^i ^ ° 


i = l 

i=l 

i=l 



^ i: 


:U; 




I" 


i 

■h 




and rearranging terms yield 


. . . •••la 

-^(pH) + V*PqH+ 2j|a^ +V *pJ? hM = -^+V. (T. 

j=l 


q) 


NS 


i = l 


(2.67) 


where 


H = h + q /2 


= hj + (qj) /2 


( 2 . 68 ) 

(2.69) 


Expanding the first two terms of Eq. (2.67) 


at ■ at 

and rearranging terms yield 


+ ^ ^ + p q .VH + H V -pq , 


+ q • VH 


) + h( 




at 


+ V • p q 


Recalling the gas continuity Eq. (2.19) 


H +Y-fiq = 0 . 


and applying the definition of the substantial derivative to the first term in 
the above expression 


aH . - „„ DH 

+ q • VH = 


at 


Dt ’ 


the first two terms of Eq. (2.67) reduce to 


Q 7~\w 1 aTIJ 

-(pH) +V.pqH = P^ . 


at 


(2.70) 
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The third term of Eq. (2.67) may be manipulated in a like manner. 
Expanding 

^ ^ ^ + pj qj . VH^ + hV • pj qj) , 

j = l 

and rearranging terms yield 

^ p (# + ? • Vrf)] + £ [rf (-^ + V . pi ?)] . 

•j-1 j-1 


Recalling the particle continuity Eq. (2.26) 


at 


+ V • ^ = 0 j = 1, NP , " J 


and applying the definition of the substantial derivative to the first term in 
the above expression 


^ +^. vni = 2 Ssi 


at 


Dt 


the third term of Eq. (2.67) reduces to 
NP, 




_L_ r . . . * .1 IN ir' . 

(pj hJ) + V . pj qj Rj] = ^ P^ 
1 ^ 1 ^ 


DH-’ 
Dt • 


Substituting Eqs. (2.70) anri f^^ into Eq. (2.67) results in 
j=l i=l 


Q 


Recall that 


rJ = hj + 
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Therefore, 


DH^ Dh^ . 1 


Dt 2 Dt 


DhJ 1 D ,-^j -"j, 

dT ■'■ZDt > 




i ( 2 ? 


Dt / 


finally, 


DHf . Dh^ .2ai 

Dt “ Dt ^ ^ Dt 


(2.73) 


Substituting Eq. (2.73) into Eq. (2.72) results in 

NP / . NS 

= If +V.(?-5+f;?i.Pi£;+V.Q (2.74) 


Equation (2.74) is the general global energy equation written in terms of the 
gaseous and particle species. 

To determine the effect of particle/gas conduction and radiation, the 

th 

particle energy baleuace for the j particle is written as follows 

^ = - kJ j4jr (rj) j (T^ - T) (r^l (T^) - j ; (2.75) 


where 


= particle heat transfer (2.76) 

2r"^ film coefficient. 


Defining the Nusselt number parameter 


qJ = H— 

<«“>Stokes 


(2.771 


where 


^^^Stokes “ ^ ’ 


(2.78] 


2-25 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 






reproducibility of the 
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the particle heat transfer film coefficient may be rewritten as 


= 


kG 


J 


Furthermore the heat quantity of the particle is defined as 

3 

= m^ ^ (r^) h^ . 


Differentiating Eq. (2.80) 


rearranging terms 


DQ-' 

Dt 

Dhj 

Dt 


. 3 

= m-* t ,r (r J) 


Ji 


Dh^ 
Dt ' 


DQ-^ 


. .3 Dt ’ 

4jr m'^ (r'^) 


and substituting the above results, with the expression for back into Eq 
the particle energy balance equation becomes 


^ (rj) j (T^ - T) - o4,r (r j) (Tj) -a^ T^j 

Solving for Dh'^/Dt in the above equation 


3A^ ' Dt 

r“ 


Dh^ 

Dt 


= 3 Uj(Tj) - aj 


(rJ)‘ 


and recalling Eq. (2.46) 


m' 


= 5 


m'^ r^ 


.4 


2 . .2 

J J 
m-' r-* 


the equation for Dh^Dt may be written as 


DhJ 

Dt 


= _ ^ gIA- (T^ - T) - . [^^(t 3) - T^l 

7j m'^ r'^ I ■* 
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(2.79) 


(2.80) 


. (2.75) 


(2.81) 


(2.46) 


(2.82) 


Equation (2.8Z) may be simplified by letting 


gJ c 


CJ = -^-2- = 


k 


Pr pfJ 


(2.83i 


Therefore, 


DhJ 

Dt 


r . .4 .1 

I C j (T^ - T) (T^) - «•’ t'^ . (2.84) 

mJ I J 


This is the geaeral particle energy balance equation. The temperature of the 
particle depends upon its state, where the state may be a liquid, a liquid in the 
process of solidifying, or a solid. The temperature is uniquely related to the 
particle enthalpy by the equation of state T^ = f(h^), tabulated. 


Substituting Eqs. (2.47) and (2.84) into the general global energy Eq. (2.74) 


yields 


^ Dt 


NP I 

|- J (T^ - T) - (T^) - t"^] + q ^ • p^ Aq^ 

I r^ I ■* 

3 = 1/ 1 


NS 

= +V.(?.q)+^q;.p.£j+V 

i = l 


Q . 


This is the general global energy equation with radiation effects. 
Expanding the first term of the above equation 


« dh _ ^ /an . -- _„\ 

^ Dt ^ \at ^ ’ 


and applying the definition for total specific enthalpy 

z ; 

the first term becomes 


DT-T 1 -A 

^ P + P ‘l * + P ^ ^ • ‘l) • 


(2.85) 


(2.68) 


at 
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Recalling the momentum balance relationship 

NP 

Y^(q*q) = q X (vxq) - i Aqf - ^ , (2.52) 

j=l 

and substituting it into the above equation yields 

NP 

p ^ = p + p q“> ^yh - ^) + pJ ,[^x (7 xq)] - . (2.86) 

j=l 

This expression may be further simplified by noting that for an arbitrary 
vector A the following identity exists, , 

A • [a X (7 X A)] = 0. 


Applying the above identity to Eq. (2.86) 

_ DH _ 8H . _ VP 


^ NP . . 

Dt = P^^+Pq.(Vh--^)-q.X:P''^'^^' ’ 

j=l 


(2.86s 


and substituting the result back into Eq. (2.85) yields 

NP NP 

1- f p^ A^ C'^ (T?- T) 


p-f +Pq .(Vh - ~)-q.^ pjAJ Aqj+J^ 

j=l j=l 


j - r p'^ |^ ^(tV - O? T^j + q^ • A^ Aq^ 


m 


AID -Si* 

= + V . (T . q) 


NS 


qj^ . Pj^T + V . Q 


i=l 


For the case, in which the flow may be described as steady state, 
adiabatic, inviscid and no body forces present the above equation becomes 

NP . . . NP 

\7H. T— > 1.1 -»-i 

pq . 


JNi-' . . . IMP 

(Vh - ^) - q , pj A^ A? + - I aJ cj (T^ - T) 

j=l j=l I 


m 


ftr (f 


. 4 


) - a j t'^ I + ? • p^ A^ A? 


= 0 (2.86 b) 
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This equation can be reduced to a simpler form by letting 


= q . Aqj - qj • Aq^ + 1 (T^ - T) + — ^ (T^) - T"*! (2.86c) 

^ rY-»J L J 


m*' r*' 


Therefore, the global energy equation becomes 


_ NP . . , 

pq • (yh - -^) - V p-’ = 0 

^ j=l 


(2.87) 


To expand the term in parenthesis in the above equation, it is necessary to 
examine the fundamental thermodynamic relationship for a system in which 
chemical reactions are occurring. 


2.3 GASEOUS THERMODYNAMIC RELATIONS 

Assumption 4 given in Section 2.1 states that the gas obeys the perfect 
gas law and is in chemical equilibrium, nonequilibrium or is chemically 
frozen. A chemically frozen gas is one in which the gas stops reacting at 
a given species concentration so that the molecular weight along a stream- 
line remains at a fixed va.lue and the ratio of specific heats is a function of 
temperature only. However, for the chemically reacting case the local 
chemical species concentrations will change in accordance with type of 
chemistry assumption considered for a respective analysis. 


In the high-temperature low-velocity regions of the flow field, an equi- 
librium chemistry calculation for the gas phase is a good assumption. Local 
residence times of the flow are sufficiently long for all chemical reactions to 
proceed to completion. However, as the flow accelerates through the nozzle 
and exhaust plume the local flow residence time becomes less than that re- 
quired for the chemical reactions to reach completion. Significant deviations 
from the chemical equilibrium conditions occur and ultimately the flowfield 
chemistry usually approaches a frozen condition. This calculation is treated 
best by including the kinetics in the gasdynamic calculation, thereby avoiding 
the problem of choosing either an equilibrium or frozen chemical analysis. 
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The development of the gas thermodynamic relations and the corre- 
sponding contribution to the gasdynamic relations will consider the chemical 
kinetics. Appendix C addresses the equilibrium chemistry analysis for gas- 
particle flows. Appendix E discusses the nonequiUbrium chemistry analysis. 
Finally, a summary of the applicable equations for chemical nonequilibrium 
and chemical equilibrium flow is presented in Section 3.2.3. 


Consider only the gas phase portion of the flow system in which there 
are NS-NP different species present. If the mass fraction of each species i 
is constant, the specific enthalpy, h, of the gas depends only on specific entropy, 
S, and pressure, P. However, for a variable composition 


h = h (S.P,Xj,X2, ....X 


^NS-NP 


and thus the total differential of h is 


dh = 


ah 

as 


dS + 


‘P,X. 


ap 


NS-NP 
dP + 


s.x, 




dX. 

. , s, p,x. " 

1=1 J 


In this expression the subscript implies that the mass fractions of all 
species are constant during the variation in question. On the other hand, 
the last term in the equation is a sum of partials in each of which the S and 
P are constant, together with all but one of the mass fractions. 


For constant mass fractions the total differential of h may be written 


as 


dh = T dS + idP . 


Therefore, 


8h 

as 


P.X; 


„ , 9h 

= T and 


S,X, 


_1 

P 


By defining the chemical potential as 

= hj^ = specific enthalpy of species i 


U- 

ax. 


S,P,Xj 
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the fundamental thermodynamic relationship for a system in which chemical 
reactions are occurring may be written in the form 
- NS-NP 

dh = TdS+ ^ + ^ Mi (2.88) 

i=l 


The above equation may be rewritten as 


NS-NP 


VP 

Vh - = TVS + 




i=l 


( 2 . 89 ) 


Since the pressure is a function of the gas state variables (p, S) 


P = P(P,S) , 


( 2 . 90 ) 


the expression for VP may be written 


VP 




Noting that 


2 _ 


■ («. 


( 2 . 91 ) 


and, 


Bridgeman’s Equation 


( 2 . 92 ) 


where a, c and R are local thermodynamic equiUbrium properties, VP 
P 

becomes 

VP = a^ vp + • ^ p VS . 


Rearranging terms in the above expression 

/c_-R\ 


VS = 


(VP - a‘* VP) , 


( 2 . 93 ) 
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and utilizing the equation of state for an ideal gas 


P = pRT , 


Eq. (2.93) becomes: 


= 

\pRT ) 


VS 


(VP - a VP) 


Multiplying through by T 


/c - R\ 


TV S 


(VP - a"” VP) , 


and rearranging terms yields 


Substituting Eq. (2.97) into Eq. (2.89) results in 

2 


Vh - 


VP 


/G V / ^ V 

i-1 


Since particle species are not considered in the chemical reactions, the 
summation term upper limit may be defined as 

NG = NS - NP 

where NG =: number of gas species present. 

Equation (2.98) then becomes 
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(2.94) 


(2.95) 


(2.96) 


(2.97) 


(2.98) 


(2.99) 


( 2 . 100 ) 


Substituting Eq, (2.100) into the global energy Eq. (2.87) 


(^- 1 


±LiL£\ 

\R \ 

Kp 

p I 


’)+ ^MiVXi - A^bJ = 0, 

i=l J 3=1. 


and expanding the dot product results in 

NG 

S 

i = l 


/C -w 2 ^ NP . . . 

(-^ - l) (q • VP - q • VP) + P 52 ^i? * ^'’^i ■ 2 ° (2.101) 

= ’ j=l 


Equation (2.101) is the expanded form of the global energy equation for flow 
described as steady state, adiabatic, inviscid and no body forces present. 


Dividing through by (c /R - 1) yields 

P 

_ o NG DX. NP 

q. VP - a-q. VP+ /g - 1 ^ /R - , Ep-' 

\ P i=l P j=l 

Recalling the species continuity equation 

DX. 


= 0 . ( 2 . 102 ) 


and letting. 


and 


dT = "^i i = l.NS 


pi bJ 

®i - ^TrTT 


NG 


'•'l = c /R - 1 • 

P i = l 

the final form of Eq. (2.102) becomes 


NP 

q . VP - a q • Vp + b| = 0 . 

j=l 


(2.25) 


(2.103) 


(2.104) 


(2.105) 
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2.4 SUMMARY OF THE GOVERNING EQUATIONS FOR STEADY, ADIABATIC, 
NONEQUILIBRIUM FLOWS OF REACTING GAS-PARTICLE MIXTURES 
WITHOUT TRANSPORT OR BODY FORCE EFFECTS 


The system of basic governing equations for any reacting flow field have 
now been derived. When the reactions between components of the gas mixture 
are known and the boundary conditions adequately specified, one should be able 
to solve the non equilibrium flow system. 

The pertinent equations derived in this section are summarized as 
follows: 

Continuity equation for steady state flow 


V*Pq = 0 

gas continuity equation 

(2.19) 

V • P^ q^ = 0 j = 1,NP 

particle continuity equation 

(2.26) 

p q. VX. - w. = 0 i = 1,NS 

species continuity equation 

(2.25) 


Momentum equation for steady state, Inviscid flow with no body forces 

'NP ■■ 

p q . vq -h VP = 0 global (gas -f particle) momentum 

equation 


where 






(2.50) 

(2.46) 




Vq^ - V 


0 j = 1,NP particle momentum equation (2.47) 


Energy equation for steady state, inviscid flow, with no body forces and no heat 
loss from the gas -particle system 

VP - a^q . vp - V p^ B^+ ijj, = 0 global energy equation (2.105) 

- . :v 
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where 


pi 1 


q . Aq^ - "qj. Aq^ + | (T^ - T) + . (Tj) -a^ t'^1 I (2 . 103) 

^ aJ m^rJ I J j 


AqJ = q - qJ 


and 


cj = ^ 

V f ^ 


J 


? . Vhj + I aJ d (Tj - T) + -4^ [ (T^)^ - 4 T^l = 

^ mJ I * 


(2.43) 

(2.83) 


0 j = 1,NP (2.84) 

particle energy- 
balance equation 
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Section 3 

THE METHOD OF CHARACTERISTICS SOLUTION 

The set of equations summarized in Section 2.4 are applicable for flow 
calculations in both the subsonic and supersonic regimes. Examination of the 
equations reveals that a closed form solution is not possible, thus necessitating 
a numerical solution to obtain the flowfield properties. The choice of the 
numerical solution is governed by the application and the flow regime of 
interest. . 

The current study addresses the nozzle-plume flowfield where the flow 
is everywhere supersonic except for those regions in the immediate vicinity 
downstream of a normal shock. Treatment of these regions are considered 
beyond the scope of this study, thus reducing the problem to one of analyzing 
a supersonic flowfield. For supersonic flow applications the method-of- 
characteristics provides a reliable method fox numerically solving the 
set of governing equations. The method is characterized by rapid convergence 
of the numerical solution and has been shown to be unconditionally stable. This 
method is the one chosen for the current study. 

A number of investigators (Refs. 1, 3 and 9) have developed solutions 
employing the method-of-characteristics. In ’'characteristic form" the 
relations are transformed to a set of differential equations which apply along 
the characteristic directions, Prozan (Ref. 9 ) developed a nozzle-plume 
solution which includes gas thermochemistry considerations . Kliegel (Ref . 3) 
and Hoffman (Ref. 1) extended the method-of-characteristics to include the 
treatment of gas -particle flows , The present study followed the work of 
Prozan in the treatment of the gas equilibrium thermochemistry considera- 
tions and that of Thoenes, etc., (Ref. 13) for chemical kinetics; These con- 
siderations have been combined with the coupled gas-particle solution as 
formulated by Kliegel. 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 


The original formulation by Prozan utilized the enthalpy-entropy- 
velocity form of the compatibility equations, i.e., independent variables are 
the velocity, total enthalpy and entropy). Prozan showed that the thermo- 
chemistry calculations could be uncoupled from the gasdynamic solution. 
Consequently the thermodynamics were calculated in terms of the independent 
variables, V, H, S and retrieved via interpolation from tables as needed by the 
characteristic solution. The routines became an integral part of the stream- 
line normal code developed by Ruo (Ref. 10) from which the RAMP code evolved. 
Consequently the V, H, S form is used in the present analysis for chemical 
equilibrium calculations. However, when kinetics or transition between the 
continuum or non-continuum flow is to be treated, use of the pres sure -density- 
velocity form of compatibility equations is more appropriate. Therefore 
development of both forms is given in the following sections . 

3.1 DEVELOPMENT OF THE CHARACTERISTIC CURVES 

The governing flow equations summarized in Section 2,4 may be written 
in the following expanded two-dimensional (or axisymmetric) form; 

Continuity equation for steady state, particle species not considered in the 
chemical reactions 


(2.19) 

pu 

^ X 

+ p V + a p + V p + - 0 

V X Y 

(3.1) 

(2.26) 

P'^u^ + 

^ X ^ Y 

+ u3pj + V V + = 0 j - 1, NP 

X Y y ■ 

(3.2) 

(2.25) 

pu (X-) +pv (X.) - W. = 0 i = 1,NG 

X y 

(3.3) 


Momentum equation for steady state, inylscid flow, with no body forces 

NP 

(2.50) + ^ (n- u^) = 0 (3.4) 

3=1 
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0 


(3.5) 


p u + p vv^ ^ ^ “ 


j = l 


(2.47) 


P'^u'^uA +p^v^u^ - p^A'^(u-uj) =0 j = 1,NP 
X y 


(3.6) 


piu^v^ ^.pjyjyj -pjA'^ (v-v^) = 0 j = 1,NP 
X y 


(3.7) 


Energy equation for steady state, Lnviscld flow, with no body forces and no 
heat loss from the gas -particle system 


(2.105) 

(2.84) 


NP . 


u P^ + V Py - a u p^ - a V p^ - ^ p^ A^B-J + 1 };^ = 0 

j=l 


(3.8) 


pjJh^ +pjv^hj. +|pjA^C^(Tf-T) - a-’ T^l = 0 j = l,NP (3.9) 

y m? ^ J 


The corresponding coordinate system used in this analysis is presented 
in Fig. 3-1. 



Fig, 3-1 - Coordinate System Used in the Development of the Characteristic 
and Compatibility Equations 
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To minimize the amount of "bookkeeping" that would be required in the 
development of the characteristic equations the governing flow Eqs. (3.1) 
through (3.9) will be written in the following general form 


L -a s — + b - 5 — + c for n = 1,N and hi, = 1,N (3.10) 

m mn ox mn 9y m 


where: 


N = 4 + 4NP +NG 


^ represents the dependent variables u, v, P, p, u^, v^, p^, h^ and X. 


a and b are coefficients, 
mn mn 


If we let 


1 = u. 


*2 =V 


*3=P 


=P 




(3.11) 


^4j+3 . ^ ~ i.NP 

</>4j+4=^^ j = l,NP 


4+4NP+i = ^i i = 


equation ( 3 . 10 ) can be rewritten as 


L = a , u + b , u + a ^ V + b ^ v + P,, 

m ml X ml y m 2 x m 2 y m3 x m3 y 


1 1 1 1 
+ a /O +b .0 +a -u +b -u +a , v +b > v 
m4 ”x m4 m5 x m5 y m6 x mo y 


+ a _ -f b - p ^ + a h b q + . • . 
m 7 ^x m 7 ^y m 8 x m 8 y 


(3,12) 


+ V. 4 j+l "i + ‘=m. 4j+l "y + • • • + *m, 4j44 ''i + '>m, 4j+4 >^y ^ ' 


^ra, 4 +4NP.+i ^"^i-x *^m, 4 +4NP + i ^i^ + form= l, N 

j y * 1x1 

m=:l 
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We then have N independent, linear nonhomogeneous equations written in N 
unknown derivatives. The N linear equations for may be combined to 
form a single differential operator by employing arbitrary multipliers and 
summing. Thus, 


N 

L = V a L = 0 

ZmJ mm 

m=l 


(3.13) 


where = arbitrary multipliers. 


Assuming that the following relation holds. 




a O' = b a for n = 1,N and m = 1,N; 

mn m dx mn m 


(3.14) 


Eq. (3.13) can be rewritten in the form of an ordinary differential equation. 
Substitution of Eqs. (3.10) and (3.14) into Eq. (3.13) yields 


N 


L = 


n 


+ a a 


dy ^ n 


+ c ._ = 0 


^ j ^°^m ^mn dx ' '^m “mn dx 9y ‘ '"m '^m^ 
m=l 


or N 


dtj) 




m=l 


a dx + a a -q — dy + c a dx) = 0 

m mn ox m mn dy ■' m m 


and finally. 


N 


L = (cr a d^ + c cr dx) = 0 n = 1, N 
/ j m mn ^ n m m 

m=l 


(3.15) 


Equation (3.15) is the generalized compatibility equation and is valid if and 
only if Eq. (3.14) is true. 

Equation (3.14) may be rewritten in the form 


O' (a ^ - b ) = 0 n = 1, N and m = 1,N . 

m mn dx mn 


(3. 16) 
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If Eq. (3.16) has a solution other than trivial, i.e., all cr^ = 0, then the 

determinant of the coefficients of a must be zero. Therefore: 

m 



a 

mn dx mn 


0 n = 1 , N and m = 1 , N 


(3.17) 


The equation D = 0 is called the characteristic equation for the system of 
equations (3.16). On expanding the determinant, it is seen that Eq. (3.17) 
is an algebraic equation of degree N and thus has N real roots. These 
roots are called the characteristic roots. 


To analyze the determinant D, begin by substituting Eqs.(3.1), (3.4), 
(3.5) and (3.8) into Eq. (3.12), and putting the results into a matrix of the 
form: 


^■11’ *^11 

^12’^lZ 

^13’ *^13 

^14’ ^^14 

^21’ '^21 

^zv^zz 

^23’ *^23 

^24’ ^24 

^31’ ^^31 

^32’^32 

a h 

33’ 33 

PL Vi 

34’ 34 

P 

cr 

^42’ *^42 

^43’ *^43 

^44’ ^44 


This yields 


p,0 

0,p 

pu,pv 

0, 0 

0, 0 

pu, pv 

0,0 

0,0 


and 

C 1 = Spy 

y 


0,0 u, V 

1, 0 0. 0 

0 , 1 0, 0 
2 

u, V -a u, 


(3.20) 


NP . . . 

^2 ^ (u - u-^) 

j=l 


2 

a V 


(3.19) 
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(3.20) 

Cont'd 


NP 

^3 " 2 (v - v^) 

NP 

<=4 = +1 - 

j=l 


Next, substituting Eqs. (3.2), (3.6), (3,7) and (3,9) into Eq. (3,12) with j = 1 (for 
the first particle) and putting the results into a matrix of the form of (3.18) 
yields 


(3.21) 


1 „ 
P ,0 

0,P^ 

1 1 
u ,v 

0,0 

1 1 
P u ,P V 

0,0 

0, 0 

0,0 

0,0 

P^\pV^ 

0, 0 

0,0 

0, 0 

0,0 

0,0 

nl 1 

P u ,P V 


where 


5 y 

= -p^A^u-u^) 

= -p^A^(v-v^) 

^8 = ^'^<^2 


A'^C'^(T'^ - T) + -42-^ ^j(T^)^ . 

m^ r'^ I j 


(3.22) 


(3.23) 


Repeating for each particle, j, and generalizing yields 


= 

''■2 


0 

Q,p\ 

U'^, 

0, 0 


pJuJ, pJvJ 

0,0 

0,0 

0, 0 


0,0 

p^u'^, 

0,0 

0,0 

j 

0,0 

0 , 0 

0,0 

pjuj,pjyj 



j = 1,NP (3.24) 
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and 


G 

C 


Substituting Eq. (3.3) into Eq. (3.12) 



6 



4 j+l 

y 



4 j +2 

= - P^ 

(u 

-uj) 

' 4 j +3 

= - pj 

(V 

-vj) 

' 4 j 44 

= p-’ '4 




j = l.NP 


m, 4+4NP+i 


0 

if 

m 

n = 4+4NP+i 

pu 

if 

m = n 


0 

if 

m \n 


pv 

if 

m = n 



'4+4NP+1 = ■ i = >.NG; 

and rewriting in the matrix form of (3.18) /ields 


pu, pv 

0,0 


0,0 


0 , 0 , 

P' 1 , pv 
0 , 0 


0 , 0 


0 , 0 
pu.pv 


( 3 . 25 ) 


m 




(3.26) 



ST? 

IIP 

(3.27) 

II 


111 


D 


irt 

fif! 


(3.28) 


f 1 


If all of the coefficients of Eq. (3.12) were put into a single matrix of the form: 


X - 


^ir ^11 
^ 21 ’ ^21 


^N.r ‘^Ni 


®^ 12 ' *^12 


^IN’ *^1N 


^NN’ ^NN 


(3.29) 


i 


m 

mi 
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where 


N = 4 + 4NP +NG ; 
the result would be the diagonal matrix 


X = 







\ 



0 

0 . . 

. . 0 . . 

. . 0 

\ 

0 

0 

X* 

Ov,. . 

. . 0 . . 

. 0 

/ 

;7 

(/ 

0 

0 

■ 0 

xf . . 

. . 0 .. . 

. . 0 


0 

0 

0 

0 . 

. . 0 . . 

. . 0 


0 

0 

0 

0 . . 

. . x^ . . 

. . 0 


0 

0 

0 

0 . . 

. . 0 . • 

• • 0 


0 

0 

0 

0 . . 

. 1 0 . . 

. . ,NP 
^2 

0 

0 

0 

0 . . 

. . 0 . . 

. . 0 


CO 


(3.30) 


Examining matrix (3.30), from which the coefficients a and b for 
® " mn mn 

determinant D will come, it can be seen that D will have the same form as 
matrix (3.30). That is: 


D = 


where 


^1 

0 

0 . 

. . . 0 . 

. . . 0 

0 

0 

1^2 

0 . 

. . . 0 . 

. . . 0 

0 

0 

0 


. . . 0 . 

. . . 0 

0 

0 

0 

0 . 

, . . 0 . 

.... 0 

0 

0 

0 

0 . 

D'^ 

. . . 1^2 • 

. . . 0 

0 

0 

0 

0 . 

. . . 0 . 

. . . 0 

^NP 

0 

0 

0 

0 . 

. . . 0 . 

• • -“z 

0 

0 

0 

0 

0 . 

... 0 

D3 


(3.31) 


|DJ obtains its elements from 
jD^I from X 2 > ®^c. 
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For purposes of simplification, let: 


S = u 4^ - V 
dx 

Sj = 

dx 

y. = ^ 

^ dx 


(3.32) 


Substituting the coefficients from matrix (3.19) into matrix (3. 17) and using 
Eq. (3.32) results in 


D, 


py' 

pS 

0 

0 


-P 

0 

ps 

0 


0 

- 

?-l 

s 


s 

0 

0 


-a^S 


(3.33) 


Matrix (3.33) may be rewritten in the form 

' li - V 

|Dj\/|;d.(-i)i-^N.. 
rH ij ij 

^ j = l ; , . . . 


(3.34) 


where 


and 


d = elements of D, 

IJ 1 

N. . = the minor of d. . . 
ij tj 


Utilizing Eq. (3.34) and letting i = 3, Eq. (3.33) becomes 
,0 


°l| ’ 


(-1)^N32 +d33 (-1)^N33 + 




Expanding terms 


D, = -pS 


py' 

0 

s 


py* 

-p s 

ps 

0 

y ' 
s 

0 

-a^S 

+ (-1) 

ps 

0 

0 0 
0 -a^S 
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performing the necessary matrix algebra 


|dJ = -pS [-pa^ S(y')^ +pS^] -1 (-p^a S^) ; 

and combining terms yields 

l^ll = p(y')^ -S^ + a^] • 

Substituting the expressions for S and y' from Eqs. (3.32) 


expanding 

Id 


„2/ „\2[ 2/^\^ 

= P P \dx; \dx/ 


+ 2 uv g + a' 


and recombining terms yields the final form of jD^ 

v2 

. 1 7. / dv \z 1. z lav 

D 


1 


= (“ S - (to) 


+ 2 uv ^ + (a^ - v^yj 


(3.35) 


The same procedure is applied to matrix (3.24). Substituting the coefficients 
from matrix (3.24) into matrix (3.17) and using Eqs. (3.32) results in 


loil a 


pjyl -pJ 


p%^ 0 

0 




0 

0 


pJs^ 0 0 

0 0 P^S^ 


j = l.NP 


(3.36) 


Utilizing Eq. (3.34) and letting i = 2, Eq. (3.36) becomes 

,0 >0 


°2 ^ *^21^'^^ 


>0 

^21 '-‘>^”22 *%3 <- 1 '® N 23 +^4 <' 


1)‘n24 
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Expanding terms, performing the necessary- matrix algebra 

-pj 0 

|p^ 0 0 


= p^ (-1) 


0 


0 pJ 


.2,3 

= - p^sj [-(pV (sV] . 


and combining terms yield 


II .3.4 

Di = (pJ) (SJ) 


Substituting the expression for from Eq. (3.32) yields the final form of |d^ 


(3.37) 


Finally, the same procedure is applied to matrix (3.28). Substituting the 
coefficients from matrix (3.28) into matrix (3.17) and using Eq.(3.32) results 


in 


^3 = 


/3iS 0 ... 0 

0 pS ! 

6 ... 0 . . ps 


= (pS) 


NG 


Substituting the expression for S from Eq. (3.32) yields 

vNG 


14 = 


(3.38) 


Since the matrices jD^j through jD^] are the diagonal elements of matrix |D|, 
the characteristic equation for the system of Eq. (3.16) can be -written as 


D 


Dj * 


D 


NP 


* Dj = 0 


(3.39) 
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Substituting Eqs. (3.35), (3.37) and 3.38) into Eq. (3.39) 

v2 


l°l = to ■ (to) + 2 IIV ^ + (a^ - v^)j p^° (u ^-v) 




dx 

and combining terms; the final form of the characteristic equation becomes: 
2+NG t ^ , v2 




+ 2 uv ^ + (a^ - v^) 




(3 


Realizing that 0; the characteristic roots of Eq. (3,40) are determined 
by setting each of the three remaining terms to zero and solving for dy/dx. 


Setting the first term to zero 


/ to 

(“to-''] =0. 


and solving for dy/dx results in 


^ i = tane . 

ux u 


(3.41) 


Setting the fourth term to zero 


- vj)^ = 0. j = 1.NP 


and solving for dy/dx results in 


= tane^ j = l.NP . 


(3.42) 


Referring to Fig. 3-1, Eqs. (3.41) and (3.42) shows that each of the character- 
istics is inclined at an angle tangent to the gas velocity vector and particle 
velocity vectors respectively. Thus we see that the characteristics are 
identical with the gas streamlines and particle streamlines of the flow. 
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The final characteristic root is obtained by setting the second term in Eq. (3.40) 


Applying the quadratic formula 


dx 


2 uv + y4 u^ - 4 (a^ - u^) (a^ - v^) 
2 (a^ - u^) 


expanding the terms under the square root 


dx 


- 2 uv ± 2 Vu^ - a'^ + a^ + a^ u^ - u^ 
2 (a^ - u^) 


and finally moving a outside the square root yields 


ISL _ 

dx 


- uv + a^ y(u^ + v^)/a^ - 1 

“ 2-2 
a u 


(3.43) 


The first term in the square root may be simplified by utilizing the definition 
of the gas Mach number 

' M = q/a 


2.2 


(3.44) 


Substituting Eq. (3.44) into Eq. {3 A3) yields 


dx 


- XXV +a.^ - 1 

2 2 
a - u 


(3.45) 


2, /^Y 


+ 2uv ^ = 0 , 


and solving for dy/dx in the following manner: 


dividing through by 2 and combining terms under the square root 


^ = 
dx 


J 2 , 2 ^ 2, 4 

- uv + wa (u + V ) - a 

“ 2~ 2 
a - u 
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Referring to Fig. 3-1, the following relations may be written; 


and 


u = 

q COS0 

V = 

q sin0 



a = 

sm ( 

a = 

q sina 

cota 

= 


sin"^(a/q) 


Substituting Eqs. (3.46) through (3.50) into Eq. (3.45) 


dx 


2 2 2 
- q COS0 sinO q sin a cotO 

2 , . 2^ 2_v ” 

q (sin a - cos 9) 


dividing through by q 


^ - 
dx 


-cos9 sin0 + sina cosa 

2 2 
sin a - cos 9 


and multiplying the numerator and denominator by -1 yields 


dy _ cos9 sin9 + sina cosa 

dx ““ 2^ . 2 

cos 0 - sin a 

To simplify Eq. (3.51) consider the following trigonometric identities. 


and 


sin(0 + a) = sin0 cosa + cos0 sina , 

cos (9 +a) = COS0 cosa + sin0 sina , 
2 2 

cos 0 - sin a = cos(0 + a) cos(0-a) 


Multiplying Eq. (3*52) by (3.53) 


— 2 — 2 
sin(0 + a) cos(0 + a) = cos a sin0 cos0 + sin 0 sina cosa 

~ 2 2 
+ cos 0 sina cosa .+ sin a sin0 cos0 


and combining terms yield 


. .2 2 

sin(0 + a) cos(0 J* a) = (sin a + cos a) sin0 cos0 

— 2 2 
+ (sin 0 + cos 0) sina cosa ; 
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(3.55) 


■ sin(9 + a) cos(9 +a) = cos9 sin9 + sina cosa . 


By substituting Eqs. (3.54) and (3.55) into Eq. (3.51) the final form of dy/dx is 
obtained 

j sin(9 + a) cos(9 + a) • la ~r 
dy _ • ^ ' _ sin (9 4 a) 

dx “ cos (9 + a) cos (9 - ot) “ cos (9 T a) ’ 
or 

^ = tan(0 +a) . (3.56) 

Referring again to Fig, 3-1, Eq. (3.56) shows that each of the two char- 
acteristics is inclined at the Mach angle to the gas velocity vector . Thus we 
see that the characteristics are identical with the gas Mach lines of the flow. 


In summary, the characteristics of the flow have been found to be the 
gas streamlines, the particle streamlines (one for each particle species) and 
the gas Mach lines. 

3.2 DEVELOPMENT OF THE COMPATIBILITY EQUATIONS FOR GAS- 
PARTICLE FLOWS 

The relationship assumed in Eq. (3.14) 


a O’ ^ = b a 
mn m dx mn m 


for n = 1,N and m- = 1,N 


is valid along each of the characteristics. Therefore, we can now substitute 
the governing flow equations into Eq. (3.14) and knowing the conditions that 
exist along each characteristic, solve for the multipliers, o^. The corre- 
sponding compatibility equations that are valid along each characteristic are 
then derived by substituting the proper multipliers into the general compati- 
bility equation. The compatibility equations will be first derived in the 
pressure-density- velocity form. This form of the equations is used by the 
RAMP computer program when the finite rate chemistry option is utilized. 
In this form, the local temperature may be solved for directly eliminating 
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the need for a time consuming iterative solution. The compatibility equa- 
tions will then be rederived in the enthalpy-entropy -velocity form. This 
form of the equations is used by the RAMP computer program when the 
equilibrium and/or frozen chemistry option is utilized. In this form, the 
thermochemical data generated, by the TRAN72 computer program may be 
used. 

3.2.1 Pressure-Density- Velocity Form of the Compatibility Equations 

Substituting the governing flow, Eqs. (3.1) through (3.9) into Eq. (3.14) 
results in 

for n = 1 ; 
or 

^1 S + ^2 ix ■ ° 

for n = 2 : 

- aiP + p(u^- v )03 = 0 

or 

O’! = ^ j (3.58) 

for n = 3:. 
for n = 4; 

(« ^ - V) CTi - a^ (u ^ - v) 0-4 = 0 ; (3.60) 


0 j = l,NP; 

j = l.NP ; (3.61) 
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for n = 4j+l: 


dx ^4j-H ^ P dx ^ > ^4j-t2 


or 


^cr. . , + (u'^ ^ - v^) o. . _ = 0 
dx 4j-H ' dx ' 4j-t2 


bepeoduoibility op raE 
OftlOMAIi PAGE IS POOE 


for n = 4j+2: 



- pi <'4j+r+P^<“^ to '^4j+3 

= 0 j = 1,NP 


or 

''4j+l = '“i to ■ ''4j+3 

j = 1.NP ; 

(3.62) 

for n = 4j+3: 


j = l,NP; 

(3.63) 

for n = 4j+4: 


j = 1.NP 


or 

<“i to - ''i> °'4j+4 = “ 

j = IVNP ; 

(3.64) 


!/ H 


ft e 






I ; 

S' S 


and finally-, for n = 4 + 4NP + i: 

p (u|j- V) <T4+4 nP+1 = “ 1 = 

or, 

<“ te - "4+4NP+i = “ i = l.NG- (3.65) 

Now, substituting the values of and from Eqs, (3,19), (3.20), (3.24), 

(3.25) and (3.28) into the generalized compatibility equation (3.15) yields 


L = (p CTj + pu (y^) du + (pu CTg) + i^2 ^ (acXj - a ua^) dp 

NP 

+ E ,[(P^ '^4j+l p’ ''4j+2> “-4j+3> dvi + (ui c. dpi 

J=1 ^J+1 


. ng 

+ (P^ o^4j+4) S ‘^4+4NP+i '^^i +1^ ‘"l ■^'"2 Z) 

i=l 


NP 


j = l 


NP 


NP . . . NP r^J j 

^3 “ ■*■ ^4 ^‘^l ~ Z P^ ®1^ ■*■ 2 1 y'" °"4j+l 

j = l j = l j=l 


(3.66) 


i 
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- P’ (“ - <'4j+2 - P' A-' (V - V-’) + fi> 4 a^.^] 


NG 1 

- S °‘4+4NP+i 
i=l 

Equation (3.66) is the expanded form of the generalized compatibility 
equation and is used to determine the compatibility equations that are valid 
along each characteristic. 


(3.66) 

(Continued) 


To determine the compatibility equations that are applied along gas 
streamlines recall that along gas streamlines 


dy V 

dx u 


(3.41) 


Therefore, Eq. (3.57) becomes: 


or 






Equation (3.59) becomes 


or 


"2 g -3 + = 


dx u 

2 dy 3 V 3 


Along a gas streamline: 


J / 0 j = 1,NP 


(3.57) 


(3.67) 

(3.59) 

(3.68) 
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Therefore, Eq. (3.63) becomes; 


(u^ f - vi) = 0 j=l.NP 


or 


^4j^l = 0 j = l.NP. 


Similarly, Eq.(3.64) becomes 


or 


^ . J) a 

dx ' ‘^4j-t4 


°^4j+4 - ° 


= 0 j = 1,NP 
j = 1,NP . 


Recalling Eq. (3.69), Eq. (3.62) yields 

0 



or 




'^4j+3 = “ j “ 


Similarly, Eq. (3.61) yields 


or 




°"4j+2 "" ° j = 1,NP . 


Substituting Eqs. (3.67) through (3.72) into Eq. (3.66) 

2 2 
L = o^ du + pu (y °’4) ^ (-a u a^) dp 


+ XI ‘"4-f4NP+i ^i V Z ^3 ^ •*• 


'4+4NP+i 
i=l j=l 
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NP . . 


NP . . . 


j=l 

and rearranging terms yields 


j=l 


12 NP , 1 

L = dvi + pu dv + ^ dP + p^ (u - u^) + (v - v^*) | dx | a. 


i i i NG , 

+ (v - v^) aj dx + (ijJj - ^ p^ B j). dx - ^ <y4^4j^p+i dx 

■ ■ i=l 


NP 


+ |u dP - a u dp + (4 Jj - p^ A^ b|) dxj 

I j=l 1 

ng , 

+ £r4+4NP+i *^^i - ^i = 0 • 

i=r 

Since the remaining multipliers, ff, are arbitrary, their coefficients 
be zero. Therefore, the coefficient of yields 


2 

PiL. 

V 


NP 


du + pu dv + ^ dP + P'^ A'^ (u - u'^) + (v - v^)| dx = 0 


or 


Since 


j=l 


pu du + pv dv + dP + P'^ A"^ |(u - u^) ^ “ v'^)| dx = 0 


j=l 


2 2^2 
q = u + V ; 


or 


2q dq = 2u du + 2v dv 


q dq = u du + V dv 


Substituting Eq. (3.75) into Eq. (3.74), and dividing through by p yields: 

NP _ 

q dq + ^ ^ P'^ A^ [(u - uV + ^ (v - v^)| dx = 0 . 

j=l 
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0 


NP 


u dP - u dp + ijjj dx b| dx = 

j=l 


or 


vj; NP . . . 

dP - 3l dp + dx - •i p^ b| dx 

j=l 


= 0 


and finally for the coefficient of 


(3.77) 


(pu dXj^ - Wj^ dx) = 0 


or 


i = l 


pu dX - Wj^ dx = 0 i = 1,NG 


(3.78) 


Equations (3.76), (3.77) and (3.78) are the compatibility equations that 
apply along the gas streamlines. 


To determine the compatibility equations that apply along Mach lines, 
recall that along each Mach line 


^ _ 


dx 


and, 


= tan (9 + a) 


dy / V 




dx ' u ’ 


(3.56) 


Equation (3,60) may be rewritten in the form: 


(U g - V) (a^ - a^ O 4 ) = 0 


(3.60) 


Therefore, 


or 


CTf - a CT^ = 0 


O’ j - a < 7 ^ 


(3.79) 


Equation (3.57) becomes 


^ (u ^ - V) 

^1 dx ^ ^2 dx 
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2 (u dy/dx - v) 


Using Eq. (3.32) 


S = u ^ - V ; 

dx ’ 


the expression for <72 becomes 

a- 


-a Of. dy/dx 


Equation (3.58) becomes 


or 


( 7 , = 


< 7 , 


(u|j- V) or. 


2 

Oi a 04 


Recall that along a Mach line 




/ 0 j = 1,NP 


Therefore, Eq. (3.63) yields 


or 


(aJg- vJ) =0 j = 1, 


® i “ liNP . 


Applying Eq. (3.82) to Eq. (3.61) 

0 


or, 





'4j+2 


Similarly, Eq. (3.62) yields 

.0 


to - ®4j+3 = “ i 
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or, 


o^4j+3 = 0 j = 1,NP._ 


(3.84) 


Equation (3.64) yields 


or 


dx ^ ' ‘^4j+4 


°^4j+4 = 0 


= 0 j = 1,NP 
j = l,NP. 


(3.64) 

(3.85) 


Finally, Eq. (3.65) yields 


(u ^ - V) a, 


4+4NP+1 


. = 0 i = 1,NG 


or 


a 


4+4NP+i 


. = 0 i = 1,NG . 


(3.65) 

(3.86) 


Substituting Eqs. (3.79) through (3.86) into Eq. (3.66) 

2 ^ ^4 dv Pu a^ V CT. dy/dx 

L = (p a^ cr^ - — g 3 2Z.v^- • ^ ^ 


^) du + — g— dv + |- 


ua4) 


dP 



u a. ) dp + I a^ a . - — ^ 


a CT, dy/dxN]P . 


y "4 ■ S 
/ NP 

fi-E 

' j=i 


^pj (u-J) 
j = l 


dx = 0 ; 


and rearranging terms yields 

r) 

2NP . . 


L = ,|pa^-a^^g)du.ea|i 


= !(' 

+ r^a^+i 


f 


^ dz) 






NP . . .- 

+ 4^1 - p^ A’^ B j 
j=l 


dx 


®4 = 0 


II 


i^r 




1. 

'mcjHi 


Hi' 


y 


Of- 


. iJ 

li 

O 


i!t ■'» 

f; I 

y- 
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CT J!» 


Since the multiplier, is arbitrary, the coefficient must be zero. Therefore, 

p . . (u - 4 g)dP dx 


2NP . . , . . 1 

+ ^ ^ | (v - v^) - (u - u^) ^ dx + ijj j dx - ^ dx = 0 (3.87) 

j=l ■ j=l 


NP 


The coefficient of du may be rewritten in the form 

2 


= pa^(-|) . , (3.88) 


The coefficient of dP may be rewritten in the form 

2 2 2 2 
a*“ dy _ uS - a dy/dx u dy/dx - uv - a dy/dx 

" ■ S dx = S ~ S 


or 


a^ dy _ (u^ - a^) dy/dx - uv 
" ■ S dx ~ S 


(3.89) 


Recalling that 


dx 


- uv 




- 1 


2 2 
a - u 


and rearranging terms 


(a^- u^) = - uv + - 1 

(u^-a^)^ = uv + a^^M^-l 


finally, 


(u -a ) _ uv = + - 1 . 


(3.90) 
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Substituting Eq. (3.90) into Eq. (3.89) yields 


^ ^ V 

" ■ S dx ■ S 


M^-l 


Substituting Eqs. (3.88) and (3.91) into Eq. (3.87) 

- — pv dn + ^ pu dv + ^ - 1 dP + dx 


2 NP . . 


NP . 


+ P'^ j^(v - v^) dx - (u - u^) dy^ + ijjj dx - p^ b| dx 

j=l j=l 

2 

dividing through by a p/S and rearranging terms yields 
udv - vdu + ‘J - 1 S dx 

^ p V p.y 


= 0 ; 


NP 


+ P^ A'^ |(v - v^) dx - (u - u^) dy - -y b| dxj = 0 


Recall that 


u = q COS0 , 
V = q sin0 , 


and dy/dx = tan(0 + a) along each Mach line. 


Therefore, S dx may be rewritten in the form 


S dx = ( u ^ - v) dx = (q cos0 tan(0 + a) -q sin0) dx 


and upon rearranging terms 


o j j n sin{0 + a) . „ cos(0 +a) 

S dx = q dx COS0 — — ^ ^ - sm0 ‘- 


cos(9 +a) 


cos (9 +a) 
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S dx 


q dx 


cos (9+a) 




COS0 sin(0 4 0 -) - sin 0 cos 


(e+qt)j 


(3.94) 


To simplify the above equation consider the trigonometric identity: 

sin|(9 +a) + = cos6 sin(0 + a) + sin0 cos(9 + a) . (3.95) 


Therefore, 


S dx = - — — sin[(9 +a) - 0 = q dx ^^9( ^ 

I J /a 


or 


cos(0+cx) ^ ■* cos(9+CX) 

sdx = + dx . (3.96) 

cos (9 +a) 

From Eqs. (3.46) and (3.4?) above the following expressions can be written: 


and 


du = - q sin9 d9 + cos9 dq 


dv = q cos9 d9 + sin9 dq 


(3.97) 

(3.98) 


Therefore, the first term in Eq. (3.93) may be written as 
2 2 2 2 

u dv - V du = q cos 9 d9 + q dq sin0 cos9 + q sin 9 d9 - q dq sin9 cos0 


or 


and finally, 


2 2 2 

u dv - V du = q d0 (cos 9 + sin 9) 


u dv - V du = q d9 


(3.99) 


Substituting Eqs. (3.96), (3.99) and (3.50) into Eq. (3.93) yields the final form 
of the compatibility equations that apply along Mach lines 


dx 


d9 Tcota ^ + il.y+-a_sia2_\ 

P \ y paVV cos(eTa)/ 

NP 

+ [(v-vVdx - (u - J)jdy - — b| dx] = 0 

3=1 ^ 
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Finally, to determine the compatibility equations that apply along the particle 
streamlines recall that along each particle streamline 


or 


and that, 


= id = taneS 


^ / V 

dx ^ u 


(3.42) 


We may immediately rewrite Eqs. (3.79), (3.80) and (3.81) which were derived 
for Mach lines. 


Equation (3.79) 
Equation (3.80) 
Equation (3.81) 


CTi 


= 


2 

= a a. 


2 

a a 


S dx ’ 


a cr , 


0o = 


(3.101) 


(3.102) 


(3.103) 


Then applying Eq. (3.42) to Eq. (3.6l) yields 


S-43+1 ° ^ = 


(3.61) 




= 0 j = 1*NP 


or 


dx 4j+l 


°^4j+l = ° ' 


Similarly, Eq. (3.65) yields 

(“ to - 


4+4NP+i 


= 0 i = 1,NG 


(3.104) 


(3.65) 
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‘^4+4NP+i " ° 1 - l.NG 


Substituting Eqs. (3.101) through (3.105) into Eq. (3.66) 

L = ^pa^cr^ + (na 


2 

a <y 


4 ^ 


4 S dxi 


dP 


NP 


+ (u o^ - u a^ CT^) dp ^p'^ u'’ ■*■ °"4j+3 

• j = l 
2 


>5 ^ - 2 NP 

. . .-1 opv a cr^ a . 

+ pj u^ j+4 J ■*■ Y s — dx 2^ 


2 NP 
a o^ 


j=l 

NP 


p^ A^ (v - v^") dx + cr^ dx - ^ 


j=l 


NP 


[-P^ A^ (u - uj) a\{v - yj) + p 

j = l 

and rearranging terms yields: 


j = l 

4j+3 ' '*‘2 '' 4 J+ 4 ] 


dx = 0 




(“-%-3^) 


^ ■ ' dP 4 dx 


NP 


+ ^ p^ A^ f(v - v^) dx - (u - u^) dy - -^ dxl + i|j ^ dx 

j=l 


NP 


+ ^^^jp j u^ du^ - p^ A'^ (u - u'^) dx| CT 

j = l 


4j+2 


NP 


+ W u^ dv^ - p^ A'^ (v - v^) dxj 0 - 4 j ^3 

j=l 


NP 


+ I^P^ u-^ dh^ + p^ dxj cr^j ^4 = 0 

j = l 
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Since the remaining multipliers, a, are arbitrary their coefficients must be 
zero. Therefore, the coefficient of <y^ yields 

- s t) ^ p-' I - 4 ij) ^ 

2 

+ P'^ A^j^(v- v^) dx - (u -u^) dy - -^b| dxl+ v|i^ dx = 0; (3.107) 

j=l ^ ^ 


the coefficient of yields 

NP . 


or 


u'^ du'^ - (u - u'^) dxj = 0 ; 

j = l 

u^ du^ = A'^ (u -u'^) dx j = 1,NP . 


Similarly, for the coefficient of 0'4j^3 

NP 


or 


^p^ U'^ dv^ - p^ A'^ (v - v^) dxl = 0 

j=l ^ 


u^ dv^ = A'^ (v - v^) dx j = 1, NP ; 


(3.108) 


(3.109) 


and finally for the coefficient of a 


4j+4 


NP 

y*] j^P'^ u'^ dh'^ + p^ i(j| dxj = 0 

j=l 

or 

u'^ dh'^ = -i|j^dx j = l,NP. 

Recalling that 

iljj = f aJ cj (T^ - T) + - 4 ~ (T^)^ 

m-^ r-^ ^ 



(3.110) 


(3.23) 
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Eq. (3.110) becomes 


dh^ = - 


aJ (T^ - T) + (T^) - aj T^ l 


dx j = 1,NP (3.111) 


Before evaluating Eq. (3.107) it is necessary to examine Eqs. (3.57) through 
(3.60) with S = u dy/dx - v non-vanishing (S / 0) as in the case along a particle 
streamline. 


Equation (3.60): 


or 


I 


S Oj - a S cr^ = 0 


^1 “ ^ ^4 * 


(3.112) 


Equation (3.58): 


or 


a, = S O’. 


a/S 


finally, 


2 

a o. 


0*0 = 


(3.113) 


Equation (3.57): 


or 


- - _LIl dir . 

S dx ■ S dx ’ 


(3.114) 


and finally 

Equation (3.59)i 


d V ^ 

3^ - CTo + S O . 0 . 

2 dx 3 4 


(3.59) 
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Substituting Eqs, (3.112) through (3.114) into Eq. (3.59) yields 


2 2 2 
^ ^ A J \ ^ ^ A 

t = 0 


[' (^) ■ ^ '^4 ’ “ 


(3.115) 


Assuming = 0, Eq* (3.115) becomes 


f- 

I V 


= 0 . 


Recalling that 


S = u,^-v 


S = u _ V 0 ; 

.J 


Eq. (3.115) may again be rewritten in the form: 


©- 


2 ^ 2 /v 

a + u I — 


9 X ^ 

2 uv — r + V 


(u^ - a^)^^ j - 2 uv^ + (v^ - a^) = 


Solving for -~t- , 

J 


9 X Ja 2 2 . . 2 2, , 2 2/ 

2 uv Hh w4 u V - 4 (u - a ) (v - a ) 

77^2 z 

2 (u - a ) 


In the above equation, either is zero or its coefficient is zero. 
Along particle streamlines: 


dy _ yJ 
dx “ j 


(3.42) 
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upon rearranging terms, 








^^,2 2 2 2-2 2 ^22 
uv 4- *wu V ~u V 4-a v +a u 


2 2 
u - a 




u V Hh y 

2 , 2 \ 2, 4 

a (u + V ) - a 

J 


2 2 
u - a 

v^ 

^ 2 
uv + a 

J, Z ^ 2. / 2 . 

W(u + V ) / a - 1 

uj " 


2 2 
u - a 

yj 

2 

-uv + a 

y(u^+v^)/a^ - 1 

uj ■ 


2 2 
a - u 


This equation has the same solution as Eq, (3.43), that is: 

= tan(9 + a) . 

However, for a particle streamline, 

i = tansj ; 

therefore the assumption that the coefficient of in Eq. (3.115) vanishes 
does not apply along particle streamlines, hence the multiplier must be 
equal to zero and Eq. (3.10 7) is not valid along particle streamlines. In 
summary, the compatibility equations that apply along particle streamlines 
are: 

i = tanej j = 1,NP 
dx 

u'^ du'^ = A'^ (u - u'^) dx j = 1,NP (3 

u'^ dv^ = (v - v^) dx j = 1,NP (3 
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and 


dh^ = - 


I cj (xj - T) + X'^jj 


dx j = 1,NP (3-111) 


3.2.2 Enthalpy-Entropy-Velocity Form of the Compatibility Equations 


Recall the pressure-density-velocity form of the compatibility equations 
that apply along gas streamlines 

NP . , . 

qdq+'~+ji^ |(u - u^) + ^ (v - v^)j dx = 0 (3.76) 


j = l 


dP - dp + — ^ dx - i £ pS aJ b{ dx = 0 

3=1 


(3.77) 


pu dXj^ - dx = 0 i = 1, 


NG 


(3.78) 


Equation (3.76) may be written in the enthalpy -entropy form using the thermo- 
dynamic relations 


X dS = dh 


dP 


NG 


H = h + 


i=l 
2 


( 2 . 88 ) 


(2.68) 


and the definition of the local speed of sound. 


From Ref. 2 it is shown that since the particle velocity and temperature 
do not change through a discontinuity and since there is a finite relaxation 
time associated with changes in the particle properties, then an infinitesimal 
weak disturbance must travel at the gas -sonic speed. Xhat is to say, the local 
speed of sound of the gas-particle mixture is unaffected by the particle and is 
identical with the speed of sound of the gas; hence, it can be expressed by the 
following 


a^* = yRX . 
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ii 


'U-i 
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II 


if « 


Hi 

u 




■ L-- ^ 


From Eq. (2.68) we may write 


dH = dh + q dq 
or 

dh = dH - q dq 


Rearranging Eq. (2.88) 

NG 

^ = dh - Tds - ^ Mi dA-. . 

i = l 

substituting for dh from Eq.(3,116) 


dP 

P 


dH - q dq - TdS 


NG 


. dX. 

1 1 


i = l 


(3.116) 


and rearranging terms yields 

dP ^ 

^+qdq = dH - TdS - ^ M^dX. . (3.117) 

i = l 

Substituting this result back into Eq. (3.76) yields the final form of the com- 
patibility equation 

NG . 

dH - TdS - ^ ^i <^i + p S J (v - vM dx = 0 (3.118) 

i=l ' ^ j=l 


Equation (3.118) may be rewritten using the definition of the local speed of 
sound. 

Solving for temperature 


T 


2 ■ 2 ^ 
q , s in q 

y R 


(3.119) 
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and substituting the result back into Eq. (3.118) yields 

0.2 NG NP . . . 

dH - dS - 2 dX. + i 2 

i=l j=l 

Next, we examine Eq. (3.78). 

Equation (3.78) may be rewritten in the form 


= 0 


or 


p q COS0 dX^ = w^ dx i = 1, NG 
w. dx 

dX. = i = 1,NG 

1 pqcosO 


By letting 


= T ■ 


Eq. (3.122) becomes 


dAT. = 


Xi dx 


i q COS0 


i = 1,NG 


or finally, 


q COS0 dX. = X. dx i = 1, 
^ 1 1 


NG 


Equation (3.77) may be written in the enthalpy -entropy form using the 
thermodynamic relations 


and 


T dS 


^1 


(^-l) 

/dP 

dp \ 

\ R / 

\P 

P / 


NG 


1 


U.. W. . 
/I 1 

c /r - 

■ 


Rearranging Eq. (2.97) 


dP - a^ dp 


p T dS 

(Cp/R - 1) ’ 
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and substituting for T from Eq. (3.119) yields 


2 2 

IT-, 2 p q sin a dS 

dP - a dp = ^ . 

^ Y (Cp - R) 


(3.1 


Replacing the first two terms of Eq. (3.77) v/ith Eq. (3.125) 

NP 


2 . 2 .. i]j 


P sin^'a dS , '^l , 1 \ ' „i , i ^ 

r(c -R) + — ^-u2-.P ^ = 

P j=i 


substituting for from Eq. (2.104) 


2 2 
p q sin a d S . R 


Y(C -R) (C -R)u 
P P i = l 


NG NP 

^i^i ” u ® 1 ~ 

j = l 


(3.1 


and using the relation =Wj^/p results in 

NG 

j = l 


P sin^g dS , 


NP 


m. 


y (Cp - R) (Cp - R) q cc 


i=l 


Dividing through by pR/(C - R), and recalling Eq. (3.119) yields 


NG 


. (C - R).^ . . . 

T dS + — ^ > jU. X. dx - P p — > pJ b\ dx = 0 

q COS0 / V 1 1 upR / j ^ 1 

i=l j=l 


(3.: 


where 


B 


J _ 


R 


1 C - R 
P 


^ • Aq^ - ^ • Aqj + I (Tj - T) + . ^ ^ (T^) - a^' t‘^J | 


(2 


r‘^ 


The final form of the compatibility equations that apply along gas streamlines 
are summarized as follows; 


Equation (3.127) 


NG 


T dS 4- 


q CO 


:os0 ^ ’z ^i^r 


dx - — P 


(C - R) 


NP 


upR 


y^pj dx = 0, 


(3, 


i=l 


3=1 
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Equation (3.118) 


dH - T dS - 


iNU INJr’ 

AV. + f (u - u'^) + ^ (v - v^)| dx 

i=l j=l I ' 


(3.129) 


Equation (3.123) 


q COS 0 dX. = X. dx i = l.NG 
11 


(3.130) 


The pressure-density-velocity form of the compatibility equations that 
apply along Mach lines are 


q^ d9 + cota ^ + (ia_sin9 + h ... qsina \ ^ 

^ ^ psi . cos (9 + a) ' 


IN X~ 

+ |^(v - v^) dx - (u -u^') dy - b| dxj = 0. (3.100) 

3- 


Substituting Eq. (3.117) 


= dH - T dS - q dq - p. dX. 

i = l ^ 


( 3 . 117 ) 


Into Eq. (3. 100) 


q^ d0 + cota (dH - T dS - q dq - ^ 


sina /6q sinO 


cos(0 +a) 


+ -4)*^+5y'p^Aj (V - vj) dx - (» - dy + , 

^ . cos(9+a)q sin a 


>{dxl 
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dividing through by q , and expanding wherever possible yields: 


-7 cota . cota ^r- -.o . cota , , cota 

d9 + — 7- dH + — T dS + dq + — 

c — a — q — c 

q q ^ q 


NG 




i = l 


— 6 sina sin9 , — 

+ dx + 


sina 4* 


1 


NP 


dx 4 


cos(9 + a) y <ipa cos(0 4a) pq 




j=l 


(v - v^) dx - (u - 


BJ dx 


q sina cos (9 4a) 


= 0 


Equation (3.131) may be rewritten in the following manner: 


Recall that 


2 2 

^ _ q sin g 




7 R 
NG 

C„/R - 1 

i=l 


and 


w. 


X. 


Therefore, the third term of Eq. (3.131) may be rewritten as 


Z 2 

cota m cota q sin a 

__ T dS = f— 

q q yR 


dS 


or 


ccta -1 no sina cosa 

^ T dS = — — =r— — dS 

^2 yR 


The seventh term may be rewritten as 


sina 


sina dx 
_ 

q pa cos(9 + a) 


1 

TcTr’^ 

p 1=1 


dx 


3 2 — 

p q sin a cos(0 4a) 
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u 


or 


sina dx 


dx 

c 7 r 


NG 


qpa cos(9+a) q sinacos(9 + a) 


(3.133) 


II 


A portion of the eighth term may be rewritten as 

(v - v^) dx - (u - u'^) dy = dx |(v - v^) - (u - u'^) 

= dx |(v - v^) - (u - u'^) tan(9 +a)J 


= dx 


(v - vj) 5l2s.(9..t«). _ slnl9,_t«), 

cos(9 +a) cos(9+a)_ 


finally, 


(v-v^) dx - (u - u'^) dy = — j (v - v^) cos(9+a) - (u - u'^) sin(9 +a)| . (3.134) 

cos (9 + a) ^ ' 

Substituting Eqs . (3.132), (3,133) and (3.134) back into Eq. (3.131 yields 
the final form of the compatibility equations that apply along Mach lines 

, cota . sina cosa dS — cota dH — 6 sin9 sina dx 
d9 + — r — dq + ■rrFi + + 


- q 


rR 

NP 


dx 


-L 2 

pq cos (9 +a) “ 






NG 


c /E - 1 £ 

^co|c.y^ -_p U1 

— 2 1 ^ 3 


y cos(9 +a) 


_ __ . _ ID' 

Hh (v “ cos (9 +a) + (u - sin (9 +a) + - 
NG 


H 


^ !' 


I i 


= 0 


i=l 


q sina cos(9 +a) 


(3.135) 


u 


3-40 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 


3.2.3 Compatibility Equations for Equilibrium and/or Frozen Chemistry 

Recall the enthalpy -entropy -velocity form of the compatibility equations, 


Along gas streamlines; 


NG . (C -R) W . . . 

T dS + dx - — VpJ dx = 0 

q COS0 Ah 'i- 1 uOR 2 —/ 1 


j = l 


and 


NG 


NP 


dH - T dS - ^ jLt. dX. + ^ + u = 0 • 

• ■ j = l 


i = l 


(3.128) 


(3.129) 


Along Mach lines; 


. cota n , sina coso; dS — cota dH -r 6 sin 0 sina dx 
d 0 + dq + -^5 — ' + -5 + — 

T q y. cos(0 +a) 


NP 

cbc V- 


— 2^- 


pq cos(9+a)j^q 


+ (v - v^) cos( 0 +a) + (u - u'*) sin( 0 +ot) + 


h 


B 


NG 


1 


q sino: 


NG 


dx 

P'" i=l 


C^l E = 0 


, cota ^ -- ^ 

±— E + 


(3.135) 


i = l 


q sina cos (0 +a) 


For any closed system which is in the state of complete equilibrium the 
terms ]S[G 


i=l 


and 


NG 

Em,x, 

i=l 

must be zero (Ref. 14) . Furthermore^ for chemically frozen flow and dX/ 
are obviously zero since there are no changes in the chemical species con- 
centration. Therefore, for the case of equilibrium and/or frozen chemistry 
the above equations reduce to the form: 
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Along gas streamlines: 


NP 


and 


- -fe- S A’ dx = 0 , 

j=i 

NP _ _ 

dH - T dS + j(u - u^) + ~ (v - v^Vj dx = 0. 


(3.128a^ 


(3,129a) 


j=l 


Along Mach lines: 


. cota j_ . sina cosa dS -7 cota dH - 7 - ^ sin 0 sina dx 
cL0 I* _ ^ ** 


- q 


yR 

NP 


dx 




J aJ 


pq cos(9+a)^ 


y cos(9 + a) 


+ (v - v^) cos (9 +06) + (u - u^) sin (9 + a ) + 


M .3 

q sma I ^ • 


135a) 


3.2.4 Summary of the Compatibility Equations for Gas-Particle Flows 

The characteristics of the flow have been found to be the gas streamlines, 
gas Mach lines and particle streamlines (one for each particle species). 

Pres sure -Density -Velocity Form (for Chemical Non-Equilibrium and Transi- 
tion Flow Applications) 

The slope of the gas streamline, 0, is given by 

= tane (3.41) 

and the compatibility equations which apply along gas streamlines are: 

NP 

q dq + ^ + i A^ f(u -u^) + ^ (v - vM dx = 0, (3.76) 

j = l ^ ^ 

tjj NP 

dP - a^ dp +-^dx - i ^p^ A'^ b| dx = 0 (3.77) 

j=l 
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and 


pu - Wj^ dx = 0 i = 1, 


NG 


(3.78) 


The slope of the Mach lines (left running characteristics and right 
running characteristics) is given by 


dx 


= tan(0 +a) 


(3.56) 


and the compatibility equations which apply along each Mach line are: 


2 ^ dP 

q d0 + cota — 


, /s q sin0 , ^1 \J— q sina \ , 

^ I — 3 .j. ^ I i .j. — a j ^ 

\ ^ pa /\ cos(0+a)/ 


+ ^^pj J(v- yj) dx - (u- J) dy - dxj = 0 (3.100) 


Equation (3.100) may be written in the more convenient form 

NP 


— , dP — 6 sin0 sina dx . 

d0 + c Ota — x- + -■>- 


dx 


2 ‘ — — 2 - . 
Pq ycos(0+a) Pq cos(0+a) 


Y'pj aj 


+ (v - v^) cos(0 +a) 


dx 


NG 


- • - 
+ (u - u^) sln(0 +a) + 


q sma 




i=l 


q sina cos(0 +a) 


= 0 . 


(3.136) 


The particle streamline direction, d\ is given by 


= ^ = tanoj j = 1,NP (3.42) 

and the compatibility equations which apply along particle streamlines are: 


u^ dJ = A^ (u-J) dx j = 1,NP 

u^' dv'^ = A'^" (v - v^') dx j = 1,NP 


(3.108) 

(3.109) 
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and 


u'^ dh'^ = - 


j A'^ C^ (T^ - T) + -4 ^ ; ■■ [(^'^ (T'^) - 4 T^^ll 

m^ r'^ I JJ 


dx j = 1,NP (3.111) 


Enthalpy-Entropy- Velocity Form (For Chenaical Equilibrium and/or Frozen 
Flow Applications) 

The slope of the gas streamline, 9, is given by 


iz _ 

dx 


tan9 


(3.41) 


and the compatibility equations which apply along gas streamlines are: 

LP 


dH 


- T dS + ^ pj A^ [(u - u'^) + ^ (v - v^) j dx = 0 

j=l 

(C - R) ^ . . . 

aJb{ dx = 0 


(3.129a) 


(3.128a) 


j = l 


where 


bJ = 


jq . A? - qj. A? 4 f (T^ - T) -f 
p' { A m-' r-' 


3cr 


(tJ) 


-4 T^ jj 


(2.103) 


A^ - -i 
^ " 2 


l/f'’ 




and 


j _ 


c-> = 




The slope of the Mach lines is given by 


^ - 
dx ~ 


tan(0 +a) 


(2.46) 


(2.83) 


(3.56) 
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and the compatibility equations which apply along each Mach line are: 


. cota 1 . sino cosa dS -r cota dH — 6 sin0 sina dx 

d0 + __ dq + -5 + 2" + 

^ ^ 'I y cos(9 +a) 


dx 

2 

pq CO 


NP r 

)s(9 + a)r^ [ 


+ (v - v^) cos(0 +a) + (u - u'^) sin(0 +a) + 


q sina 


(3.135a) 


The independent variables are x, y, q, 9, P, p, v\ h-^ and in the 
Pr es sure-Density- Velocity form, and x, y, q, 9, H, S, v^, h^ and p^ in the 
Enthalpy-Entropy- Velocity form. 6+4NP unknowns are required to corn- - 
pletely define the gas -particle flow at a given location in the flow field. The 
above equations provide only 6+3NP compatibility relations. This is because 
a compatibility relation does not exist for the particle density and results 
from the assumption that the particles do not contribute to the pressure acting 
on a control volume in the gas -particle system. As originally shown by Kliegel, 
(Ref. 3), this can be resolved by utilizing the following incompressible flow re- 
lation to obtain the particle density. 

. .5 . . .6 .) 

(yJ) dy^ - v"^ (y^) dx*^) (3.137) 


dm^ = (2ir)^ 


and 

6 takes on the values 

5 = 0 for two-dimensional flow 

1 for axisymmetric flow 

3.3 FINITE DIFFERENCE SOLUTION OF THE COMPATIBILITY EQUATIONS 

To solve the system of compatibility equations it is first necessary to 
write these equations in finite difference form. Notice that the equations de- 
scribing the particle properties were derived for a discrete particle. To 
account for more than one particle size or species, n discrete particles are 
allowed in the numerical solution. The streamline of each additional particle 
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id an additional characteristic curve; and the solution is obtained by applying 
the particle compatibility equations independently for each particle along its 
streamline. 

The equations are cast in a form suitable for the calculation of an 
interior point in the flow field. For special points, such as along the nozzle 
axis, the nozzle wall, the limiting particle streamline and a free boundary, 
certain boundary conditions in the flow field are known which allow some of 
the equations to be discarded. The characteristic net for an interior point 
is illustrated by Fig. 3-2. 

The difference equations for the gas-particle system of compatibility 
equations in coefficient form is given below. 


The slope of the gas streamline, 9, is given by 


5-3 = tan9- ^ 
Ax 5 , 3 


(3,138) 


and the compatibility equations which apply along gas streamlines are 


Pressure-Density-yelocity Form for Finite Rate Chemistry 




(3.139) 


^^5-3 ■ ^5, 3 "^5-3 ~ ° 


(3.140) 


Pc 1 'ic o AX. - w. Ax_ „ = 0 i = 1,NG ; 
5.3 5,3 ig_3 3 5-3 


(3.141) 


where the barred values are averages over the step length, and the coefficients 
are defined as follows 
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IN ir 

E H, 3 3 [<"5, 3 - “5, 3> 3 * '^5, 3 -"'s, 3’*“®5, 

Hr- O i _ 1 _ •* 


5 , 3 j = l 


- 1 NG 


T, 3*15 3 -Z 4 . 3 4 , 3 3 

^5,3 J i=l ^ j = l 5 , 3 qg 3 


AL is the step size along the streamline. 


Enthalpy-Entropy-Velocity Form for Equilibrium and/or Frozen Chemistry 


^^5-3 = J5.3^S5_3-C3 


AS5_3 = C2 


^2# 


1 al, 3 . 

2 3 4 3 z: — 


The s,F6pe of the Mach lines (left running denoted by subscript 2 and 
right running denoted by subscript 1) is given by 

^1.2 = tang, 2 C 


where 


^1, 2 "" ®1, 2 ^1,2 “l, 2 . ■ 
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and X takes on the values 


for interior solution 

for lower boundary solution 

for upper boundary solution 

for interior solution 

for lower boundary solution 

for upper boundary solution. 


-V" 


The compatibility equations^ which apply along each Mach line are: 


Pr e s sur e - Den sit y - V elocit y F orm 


^® l ,2 ■‘■^ 1.2 ^^ 1,2 ~ ° 


( 3 . 145 ) 




where the coefficients are defined as follows 

cosa 


Q 


1,2 


1^2 


' 1.2 




^X2 ‘^°=Pl.2 




and 


Cl 


1 . 2 


NP . . r _ 

23^1,2*1.2 -'^1,2 ‘^1.2’ •=°=Pl,2 + <"l,2-“l,2> *‘"l>!,2 

j -1 I 


"‘ 1.2 

1l, 2 ^‘"“l,2l 


ZilX 


Pl,2‘»l, 2 2 2 =‘""l, 2 '““P 1, 2 Si 


LA 


NG _ 

E ]I X: 

M 2 ^1 

i=l 


Ax 


LA 


^ - 1 


R 


1.2 
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Enthalpy-Entropy-Velocity Form “ 


2 ±«I, 2 ^1l. 2 1, 2 + «l, 2 2 °1. 2 i <=■ 


where the coefficients are defined as follows 


"^““1,2 ^1, 2 


sinOj 2 cosOj 2 A^j 2 

2 ^1, 2 


12 “ — —2 

slnaj.^q, 2 


^^‘“^r.2°‘”°l,2^1..2 

n,2‘=?=Pl,2 


2 M, 2 ± <'' 1 , 2 "i, 2> 2 + <“1, 2 -“1 

I j=> t 

I - i 1 - -T‘-^ - ' 

■ «1,2 *‘“l,2 |Pl,2«f,2'°‘'0l,2 


Finally , the particle streamline direction, 9^, is given by 


4-3 = tanej^ j = 1,NP 

i/J ^ f J 
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o 


and the compatibility equations which apply along particle streamlines are: 


^“4-3 

11 

'^1.3 

j = l.NP 

(3.148) 

■'^''4-3 

= 3 

^4-3 

j = l.NP 

(3.149) 

^4-3 

= ^4,3 

^4-3 

3 = l.NP 

(3.150) 


where the coefficients are defined as follows 


czi - 

= A^ 

<"4 

t.11 

4, 3 

4, 3 


“i: 

C3^ o 

= aJ 

("4 


4, 3 

4.3 


ul 


vi ,) 


O * * * A 

“ 1 ^ 4,3 ^ 4 , 3 <■^ 4 . 3 •'^ 4 . 3 ^ ^ 3~~J 


"™ 4 . 3 ^^ 4 . 3 


3 (T^ ■ 

4. 3 '^4.3' °4. 3 4. 


The difference form of Kliegel’s incompressible, flow relation is 

AmJ 2 = <yi,2> M, 2 ' ^1, 2 <yi, 2>* 


An iterative solution is employed to determine the properties of the flow 
-at the new point. For the first pass of the solution the barred values are ap- ^ 
proximated by the conditions at the known points; e*g-, 

01 = 


After the appropriate set of equations has been solved a new estimate of the 
barred values is made, e.g., a x o 

e. 


1 2 

■ ■ s 
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The iterative process is continued until the desired convergence is 
reached . • ■ . .. - *. 

The mechanics- of the numerical solution for a typical point are quite 
involved. Section 6 presents the mechanics of the numerical solution for the 

different types of points encountered in a typical gas -particle flow problem, 

3,4 DEVELOPMENT OF THE PARTICLE DENSITY RELATIONSHIPS 
Kliegel's incompressible flow relation - - 

dmj = (2 jt)^ 

will now be used to develop the particle dehsity relationships required to 
solve for the particle density at each of the different types of points encountered 
in a typical gas -particle flow problem. The p^jticle density relationships will 
be derived for an interior point, a lower ‘boundary point, a particle limiting 
streamline point and an upper boundary point. The particle density relation- 
ships for the remaining type points (such as a shock point) will not be developed 
because each is subject to the same boundary conditions as one of the above 
four cases and therefore yield the same results. 



3.4.1 Interior Point Solution 

Consider the characteristic net for an interior point 


• ♦ w : 

u^ (y^) dy^ - (y^) dx'^ i 


(3.137) 


Normal 





m 


i - ^ ^®ft“ Running 

^ Characteristic Line 

^ “~r — 












m 


Right -Running 
Characteristic Line 
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Points 3 -2: 


'if ‘P2°2'^**3“3><y3y3 ->2^^ llllli' ,1 ,, ,_ 


and points 1 - 2: 


•^1,2 ~ 


gnPiUi + P2U2) (yijrj - 72^2) 


(^ri + ^2^2^ 1 + ^2^ 


(Xj.X2) 


(3.155) 


The partlGle density relationship at point 3 may now be. determined by- 
substituting Eqs . (3 . 1 53) , (3.154) and (3 .1 55) into Eq. (3 . 1 52) and Solving for p ^ . 
Equation (3,152) becomes: 

(PiUi +P 2 U 2 ) (yjyJ - Y2yl) - (PjVi + P 2 V 2 ) {y\ + y^) (^j - X 2 ) 

= (Pl"l +P3«3> (yiYi - yjy^) - (PjV^ + P3V3) (yj + y|) (Xj - X3) 


+ (P2U2 + P3U3) (y3y3 - y^y^), 


^ 2^2 - ^ 3 '^ 3 ^ 2 ^3 ■ * 2 ^ * ^ 
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To solve for Py 
expand Eq. (3 .156) 


PlUi (yijr^ - 7272) * Pi^ <^1 + (*i * ^2^ + Pz'^Z ^^1^1" ^^ 2 ^^ 


- P2^2 (7i + 72> (^1 -^2^ = P3 U3 (7i 7i ' 73 73) 


- PlVi (JJi + 73 ) (^1 - ^ 3 ) + P3"3 (y3i^3 * ^l"l <^1^1 " ^3^ 

+ P 2 U 2 <7373 “ 5^24) - P^Z <^2 ^ <^3 ■ ?"2^ ■ Ps^3 <^1 + 


P3'^3 ^^2 ^3^ ^^3 ■ *2^ ’ 


combine terms to form coefficients of Pj,P2 and Pj 

Pll«l(7i7i - 7272) - '^i <7i + 7^) (?i ‘ ^2^ + ^2 1"2^^1^1 ' ^2>^2^ 

■'.■'■■'■ ■ ■■.'•. ■ ■■'.■’ ■ ■ ■:■■ . ■' ■ ..'. '■ . ' ■ ■■ ■' - ■■. 

> 6 , 5\ / \1 ^ \ / 6 6 * 5 6y 

- V2 (7i + 72) <^1 - ^2)] = ^3 

- '^3<7i + 73) (Xj - x^) - v^(y2 + 73) (^3 - X2)] 




- V2(72 + 73) (X3 - X2)] ; 


isolate p.~ and its coefficient 


P3 «3(7 i 7 i - 7272T - V3[lyj^^ (x^ - X3) + (y2 + 73) (X3 - X2)] 


j T 73/ \-s-j “ -»-3; T v >2 ^ iy v-»-3 -^2' 


if - r , 6 6, , 6 ^ 6. , , , 6 • 6, 

/ =Pl [«i(7i7i - 72^2> - ^^l^^l <^1 ■ *2t- "l^^l^ir ^3^3> 
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+ v^(y\ + yI) (xj - X 3 )] +p^2(yiyi - yzyl^ " 




combine terms on the right hand side of the equal sign 

= ^i|"i<y3i^3 Vyzy2) - - ^s)]} 


+ P 2 k(yi^i - ^35^3) - ^2 


and finally solving for 


ki + y^) (^1 - ( 5^2 


^^3 ■ , 6 6 ' 1 / 5 . 6 w i . / 5 . / .l [^lrry 3^3 ‘y’ 2 ^ 2 ' 

, '^S^yi^l-yZ^z) ■ '^ 3 ri‘*'^ 3 ^ (Xj-X3) + (72 + 73) (X3-X2)] I 1 

> - Vj[{ 7 [ + yk<^l “V ' " 

-V2|(yi +yk<^l ■ '‘zV- + <^3 ( 34 ^ 

3 .4 .2 Lower Bou^ Point Solution 

The lower wall point and free boundary point are the two types of lower 
boundary points encountered in a typical gas -particle flow problem. To derive 
a relationship for the particle density at a lower boundary point it will be 
assumed that: (1) the lower bounda.ry runs parallel to the x-axis (for two -phase 
only), and (2) the particle streamlines run parallel to the lower boundary when 
in the vicinity of the lower boundary. These assumptions are valid in nearly 
all of the flow problems encountered. 
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Consider the characteristic net for a point on a lower boundary 


Normal 



Right-* Runniifg % 
Characteristic Line 




Lower Boundary 


Applying the principle of conservation of mass to the V^ve system, we may 


“ 1-5 = ”1-3 +” 3-5 • 


(3.f58) 


Furthermore, from the above assumptions m^ v| and v| may be set equal 


to zero. 


Equation (3.1 5 8) then becomes 


~ -1-3 


(3.159) 


Each term in the above equation may be solved for by applying Eq. (3.137)'be- 

'.O' ■ . ' . ' . 

tween the proper points. 


Point s 1- 5: 


JiPjUj +P5U5) 


- Piv, yi(*i -*5) 


(3.160) 


Points 1-3; 


5 |^^l"l 




- ^I'^l >^1 (^1 '^ 3 ) 


(3.161) 
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The particle density relationship at point 3 may now be determined by- 
equating Eqs. (3.160) and (3461) and solving for pj, 

Equation (3.159) becomes; 

(Plbi +P 5 U 5 ) y^yj - ZPjViyi (X| -Xg) = +P3"3) ' ^PiVi/J (Xj -X 3 ) (3.162) 


To solve for 'Pj; " 

■ ' ■ ■ ■ • ' . \ ° - 

s\ ' ' • ’ ■ ' . ' ^ ■ 

combine tbrin's of Eq. ( 3 . 162) to form coefficients of p p^ and p^ 

— ■ ■ ■: ■ ■ ’ ■ - .0 - 

<* 5 '* 3 >] -^ P 5 “ 5 = P 3“3 

' . ' ' 0 ■ 'rt- ■ ' ' - ■ ' 

and then divide through by u^y^y^ 

P5U5 - 2PiVi (X3 -X5) 

Po = - — — — ■ . 


(3.163) 


3.4.3 Particle Limiting Streamline Point Solution 

Consider the characteristic net for a point on a particle limiting stream- 
line-' ■ ' : 
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Applying the principle of conservation of mass to species j in the above 
system, we may write ** 


^^5-4 = 


3-4 


(3.164) 


By definition of a particle limiting streamline, the entire mass of 
particle species j is contained in the streamtube formed by the particle 
limiting streamline for species j. Therefore m^ ^ hiay be set equal to zero. 

Note that this condition applys only to particle species j; other particle species 
may pass through the particle limiting streamline for species j. Equation (3.164) 
then becomes 


IJ 


n t, 


Tf'. 

il i 


^3-4 ^ ^5-4 • 


(3.165) 


ii i 
I 


li f 


Each term in the above equation may be solved for by applying Eq. (3.137) be- 
tween the proper points. 


Points 3-4: 

rii3-4 = (2ir)^ 
Points 5-4; 

^5-4 " 


5 6 6 6 

(P3U, +p^u^) { y ^ Yy - Y ^ Y ^) (P3V3 +p4'^4> (73 + 74) 

J (x3-k^) 


(3.166) 


§ ? 
1. i 


(P5^5+P4«4) (7575 -7474) (P5'^5+^4^4> (^5+^4^ , 

2 .6 ■ 2 ^5 ^^5 ‘^4' 


2 


(3.l67)j 




The particle density relationship at point 3 may now be determined by equating 
Eqs. (3.166) and (3.167) and solving for p^. Equation (3.165) becomes: 


6 6 6 6 
(P3U3 + p^u^) (y3y3 - 7474) - (P3V3 + P4V4) (V3 + y^) (X3 -x^) - 

(P5^5 +P4'^4> (7575 - 7474) - (P5V5 +P4V4) (75 +74) (X5 -X4) . ' 


(3.168) 


i S 
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$ 'T-'vaj’^irvr. 


To solve for p^: 


expand Eq. (3.168) 


P3U3 {y^vl - y^vl) - P3V3 {yl + yl) (X3 - X4) +P4U4 (y^yl - y^y^) 


J 


■ ^ 4^4 <^3 ■*■ 5^4) <""3 ■ ^4) = Ps '^5 <^^5^5 " ^4^4) "'"4^ 


6 5 6 6 

+ P 4 <^4 (^5 Vs ~ ^4 ^ 4 ) ■ P 4 ^4 (Vb ■'■^4^ ^^5 "^4^ 


combine terms to form coefficients of PyP^ and p^ 

^3t"3 - ^4^^ - ^^3 ^A + ^ 4 ) <^3 - ^ 4 )] = ^b['"5 ^^5^5 ‘ ^4^4^ 

-^s^As + yl) (^^-^^i^+pX^iy^yl-y^yl) - '^ 4 [<^ 5 ’^ 4 ) 


(73 + y|) (X3 - X4) 




and finally solving for p^ 


P 3 = 


1 


'^^iy^yl - ^ 4 / 4 ) - '" 3(^3 + ^ 4 ) (^3 ■ ^ 4 ) 


^ r / 5 6 , 

^5L"5^r5^5 - ^4i^4) 


Vb iyl + 74) (^5 - X4)] +^4p4<y5^5 ■ ^3^3) - ^4[^A^ + ^4> <^5 ‘ ^4> 

(73 + 74 ) tx 3 - x^)j . (3.169) 


3.4.4 Upper Boundary Point Solution 

The upper wall point and free boundary point are the two types of upper 
boundary points encountered in a typical gas-particle flow problem. Consider 
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the characteristic net for a point on an upper boundary 



Applying the principle of conservation of mass to species j in the above 
system, we may write 

’^ 5-4 = ^^ 5-3 ■*•^ 3-4 • 

Unlike the assumption made in the lower boundary solution, the particle 
streamlines are not restricted to run parallel to the upper boundary in the 
vicinity of the boundary. Therefore, the particle mass flow rate m^_2 is not 
necessarily equal to 25 ero. Physically speaking, mg_ 2 niay be considered to 
be equivalent to particles "sticking" to the upper wall or pas sing through the 
free boundary. 


The particle density relationship at the upper boundary is derived by 
following the same mathematical procedure used in the solution for the particle 
density relationship at an interior point. The final result is 


^3(^5^5 ■ ^4^4^ - V3 [(^5 ^3^ ^^5 ’^3) + (4+^3) (^3 "^4ill 




- 


(^5 +y4) ("^5 - ^ 4 ) - (4 + 4^ (^5 ■ ^3^ 
(^5 + 4 ) ^^5 '^ 4 ^ ■ ^4 ■*■^ 3 ^ ^^3 "^ 4 ^ 
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Section 4 

THE OBLIQUE SHOCK SOLUTION 

4.1 DEVELOPMENT OF OBLIQUE SHOCK RELATIONS 

Figure 4-1 illustrates a streamtube passing through an oblique shock 
wave. This wave, which is extremely thin, will cause an almost instantan- 
eous rise in pressure and temperature. For some distance downstream of 
the shock wave (in a reacting gas) a non- equilibrium zone will exist followed 
by a return to chemical equilibrium. The following analysis discusses the 
fluid flow properties in such a way.that the non- equilibrium process need not 
be specified in order to arrive at an exact solution for the gas properties. 

Consider the control surface as shown in Fig. 4-1. Applying the princi- 
ple of conservation of mass, we may write 


Similarly, conservation of momentum gives 

2 Pi’^1 2 — / P2*^Z^^ 

-(Pl + Piqi> Aj M to +(P2 + P2q2>A2 ~ tan(. -i ) 

+ f pn.e?* - /pn.e®= 0 

N.E. N.E. 

where the subscript N.E. applies to the non- equilibrium zone. Since each 
streamline locally undergoes the same process the last two terms of the above 
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and 


19*'= sin(<- 6 )ns+ cos(f- 6 ) ts 

Ivf = -cos(e- 6 ) ns + sin(e- 6 ) ts . (4.4) 

So that, after substitution of Eqs.(4.3) and (4.4) and setting each component 
to zero, Eq, (4.2) becomes. 


2 2 
p^q 2 cos(f- 6 )A 2 - cose = 0 


(4.5) 


and. 


' 2 \ . / 2 \ 

1 sin e 4 co se ~ {^2 ^2 *^2 1^2 ® ^^(^ ” ^ 


P2A2Cos(f -6 ) 
tan(f- 6 ) 


= 0 . 


(4.6) 


But from geometry it can be seen that 


~^2 _ sin(<- 6 ) 
Aj^ ~ sine 


(4.7) 


After substitution of Eq. (4.7); Eqs. (4.1), (4.5) and (4.6) become 


P 2 q 2 sin(e- 6 ) - p^q^ sinf = 0 


(4.8) 


2 2 

cos(c- 6 ) - P^<li sinf cose = 0 


(4.9) 
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o 


z z zz 

P 2 + P 2 I 2 {(-6) - pj “ Pj<lj sin < = 0. 


(4.10) 


The above set of relations contains 6., p^, p^i q^ ^ts unknown quantities, 

but ■ , •. . 


pz — p(®2^ ^2^ ’ ^Z ^ P(®2* ^2^* 


(4.11) 


Therefore, if one variable, say < , is taken as an independent parameter the 
remaining unknowns by an iterative solution. These 

equations are, of course, formally the same as the ideal gas solution. The 
difference lies only in the variation of pressure etc., with entropy and velocity. 
It is impossible to determine the location of the new equilibrium shock point 
location without a detailed description of the non- equilibrium reaction process. 
It will be assumed therefore that this zone is thin and that no significant errors 
are introduced by letting the downstream physical location lie on the upstream 
location. 


It has been pointed out previously that the characteristic equation 
derivation was based on neglecting transport properties and as such is 
necessarily restricted to continuous regions only. The oblique shock wave 
relations derived her>j then are patching lines between the continuous regions. 


4.2 ITERATIVE SOLUTION OF THE OBLIQUE SHOCK RELATIONS 


Rearranging Eq. (4.8) yields 


sin(<f- 6) = 


Pj^q^ sine 


Pz^z 


while squaring both sides of Eq. (4.9) and substituting the above relation 
yields, after simplication 
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it 


I 


i 


s 






m 

ii 


I 


I 


m 


i 


*4^ 


I f 

I 


M 


US 


1 


m 

i. 


q 2 -q 


/^\ ^ . Z . 2 

1 \P2/ 


sin € + cos < 


_1^ 

2 


= 0 ; 


and Eq. (4.10) becomes 


2 2 
P2 + P.iqi sin c 


P2 


- 1 


- p, = 0 


In functional form Eq. (4.12) and (4.13) are simply 

/ 

G ( s 2 > q^ ) = 0 
^2 ^^2 ^ ^ 2 ^ ~ ^ “ 


Furthermore , 


9G a G 

= - 5^7 '*‘12 + “*=2 

aG^ 


The partial derivatives may be written directly from Eqs. (4.12) 'S.nd 
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(4.12) 

(4.13) 

(4.14) 

(4.15) 
(4.13) 



2 



^ =PZ-^ (^Pz) qi sin€ P 2 9 -^ (inp^) (4.16) 


8 /. /N /^1 \ ^ -2 8 /. \ 

^ =P2-5^ ('"P|) -\^j P2 8i^ (*“P2) • 

. . i 

n'-.. 


Rather than calculate the partial derivatives numerically by perturbing the 

' . 'i ^ . 

functions ilnp^ and jfoip^, approxirikte values for these derivatives will be found 
by assuming that locally the gaf?behaves ideally, that is to say 


SRj _ 8 R 2 _ 8 ^ 

9 8q.2 


_ ®Po 2 . ^ 2 . . 

^^2 * ^^2 


so that 


and 


9 {Aip 2 ) 9 (fnp 2 ) _ ^ 

^®2 ^2 


(4.17) 


^ ^^PzV ^ y 
9 ^2 ^ ^ ^2 


(fnp 2 > 


^ 2^2 

Pz 


(4.18) 
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WritiiifT Eq. (4.1v5) in finite difference form: 




and 


(n+l) 

^2 ■ ‘-'Z “ aq^ 


(n) 




Since the root’Gj = = 0 is desired, are set to zero, 

resulting in 


9G 


in) 


= s<-)+ G<">-~^-G 


(n) _i 

2 aq^ Sq^ 


9G,^^^\ /vaG,^"^ aGo^"^ aG^^”^ aG,^*^^' 


a s^ . aq^ a Sp aq- 


( 4 . 19 ) 


and 


(n+l ) (n) 

q - q 


i\ 9G, 


(n) 


M+l) _ In) 


- s' 


/dG. 


(n) 


8q, 


(4.20) 


The iterative solution using Eqs. (4.19) and (4.20) is continued until the 
desired convergence of and is reached. The solution is completed by 


r . • -1 

6=6- sin 


PiT 


1^1 


sine 


(4.21) 


4-7 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 


UmODUCIBILITY OF THE 
OmcJtNAL PAGE IS POOR 


The first guess to start the solution is an ideal gas solution to the set of 
equations. If it is indeed an ideal gas under analysis the first guess is 


6 = 


tan 


tanf 


1 

sine 






_ , cose 
I ^2 ~ cos(«-6) 


and 


R, 


°2 = =1 +CT 


fn 


2y,M^ sine - ( Vj - 1) 


+ y^fn 


tan(< -6) 
tane 


( 4 . 22 ) 
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Section 5 

EXPANSION CORNER - PRANDTL -MEYER FAN 


In some cases the ''W may be required to negotiate a sharp expansion 
turn. The problem becomes two dimensional at a sharp corner (it is impos- 
sible to conceive Of an expansion corner on an axis of symmetry) and may be 
treated with a Prandtl-Meyer expansion. 



Fig. 5-1 - Nomenclature for Mach Wave Analysis 


Since a Mach wave will support pressure changes only in a direction normal 
to itself (Ref. 11); 


dv = qd6 
du = dq 



5-1 




or 



( 5 . 1 ) 


The solution to Eq. (5.1) is a straightforward numerical integration for 
the case of a known final velocity (free -boundary case). If the turning angle 
is known, however, and the final velocity is not known, an iterative solution is 
necessary to determine the upper limit. 



(5.2) 


In the mesh construction to be discussed later a fan of rays must be gen- 
erated to allow a numerical description through a large turning angle. The turn- 
ing angle is subdivided into a number Of small turns, each of which is integrated 
numerically. Corresponding to each of these small turns is a Mach wave or 
characteristic line emanating from the corner. 


Section 6 


NUMERICAL SOLUTION (MESH POINT CONSTRUCTION) 

(i. 

The calculations described previously' are point or small region sol^r- 
tions. This section presents the mesh construction and calculation procedure 
required to describe the eiitire flow field in a typical gas-particle flow problem* 
The general principles adopted in the numerical solution are similar to those 
used in the method -of -characteristics solutions. The major difference is in 
the technique of mesh construction; here, the calculation proceeds along 
normals to the streamlines instead of along characteristic lines. 

To begin the problem of describing the entire fl^ field all necessary 
boundary conditions must, of course, be supplied. Ih^'addition, a start line 
whose properties are completely defined must be designated. 

Figure 6-1 illustrates a flow field in which there are ho discontinuities 
and in which the mesh construction is terminated when the region of interest 
has been computed. 

Figure 6-2a presents the mesh construction required to solve for the 
properties of any interior point in Fig. 6-1. The J -line is considered to be 
the known normal line. The J-line maybe the input s|artline or^he line just 
calculated. The K -line is the new normal line to be calcula.^ced. Calculation 
always starts from the lowest point and moves upward along the K -line. The 
first point on the new normal is located by extending the right-running char- 
acteristic from a properly selected point (based on desired^'mesh size) on 
the preceding known normal to intersect the lower boundary based on the 
flow properties of that point. The rest of the points including the last point 
(boundary point) of the normal are located by extending the normal to the 
local streamline to intersect the next strearriline. (A boundary is also a 

I, ■ , , 
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Nozzle Wall 



Fig. 6-1 - Basic Mesh Construction for the Streamline -Normal 
Two -Phase Numerical Solution 
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J, 





E 


if 


streamline.) Once the location of the new point is known (point 3 in Fig. 6-2a), 
the characteristic line(s) and particle streamline(s) can be made to intersect 
the preceding known normal (points 1, 2 and 8 in Fig. 6~2a). The properties 
of points 1, 2 and 8 on the known normal can be interpolated from the known 
points. The properties at point 3 can then be calculated by solving the applic- 
able compatibility equations. For an interior point, both the right-running 
and left- running characteristics are used to compute the velocity and the flow 
angle at the new point (point 3 in Fig. 6 -2a). For a boundary point, either 
the flow angle or the velocity can be found with the boundary condition; there- 
fore, only one physical characteristic is needed depending on the relative 
location of the boundary with the main flow. For shock points and Prandtl- 
Meyer expansion calculations, the mesh construction is appropriately modified 
to accommodc’te each individual case. The general method of computation, 
however, remains the same as described above. 

The mesh construction and calculation procedure foi' an interior point 
solution forms the basis for the solution of the flow properties at all other 
types of points in the flow field. For that reason, the mesh construction 
and calculation procedure for an interior point will be first described in 
detail, then, to further aid in the understanding of the numerical solution, 
presented in flow chart form. Using the knowledge gained in the discussion 
of the interior point solution, the calculation procedure for each of the re- 
maining cases will be described. Any important problems that may arise, 
or other points of inte-’:*est involved in the calculational procedure will be 
discussed. 

' y- 

6.1 INTERIOR POINT 

This section presents the details of the numerical solution for the gas 
and particle flow properties at an interior point of the flow field being analyzed. 
The solution is outlined in a step-by-step procedure and is for the case when 
particles are present at all points which enter into the calculation of the flow 
properties at the new point. 
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Consider again Fig. 6-2a and assume that the flow properties are laiown 
along the J-line. The J-line can be the input start line or the line just calcu- 
lated. The point 7K will also be a known data point since it is the last point 
calculated on the K- line. 


Before the calculation procedure can be started a decision must be 
made to use either the enthalpy-entropy-velocity form or the pressure- 
density-velocity form of the eompatibility equations to solve for the proper- 
ties at the new point. This decision is based on the type of chemistry being 
considered in the flow field. If finite rate chemistry is being considered, 
the pressure -density-velocity form of the compatibility equations is the more 
convenient form to be used. If equilibrium or frozen chemistry is being con- 
sidered, the enthalpy -entropy- velocity form of the compatibility equations is 
the more convenient form to be used. In either case, the method of calcula- 
tion is the same. For purposes of discussion, it will be assumed that the 
chemistry of the flow field being analyzed may be considered to be in chem- 
ical equilibrium. 

Step 1; For the first pass of the solution, the properties at 3K are set 
equal to those of the base point 5 J. 

Step 2: The physical location (x^, y^/ of the point 3K is found by inter- 
secting the normal from the point 7K with the projection of the streamline 
from the point SJ. 


Step 3; The intersection of the characteristic lines with the J-line is 
now found. These lines will in general fall between two known data points 
and an interpolation scheme is employed to find the flow properties at these 
points. The interpolation is of the following form: 





+ ^ <^ 6 . 5 



< ' ^ .;6-5 ■■■ : ^ : f ': 
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where 


F = interpolation factor. 


For the gas flow properties P assumes the following 


P =q,e,S,H 


T .r= f{q,S,H) 
p =f(q,S,H) 

U =f{T,H) j 

Cp = f H) J._ tabulated 


Pr = f (T,H) 


while for the particle properties P assumes 


with 


P = u^ , , h^ , p^ 




f(h^)j^ 2 tabulated 


^ ^ 1 , 2 ^ 1, 2 1,2 




fj 2 = 2 

2 = f (Re^i 2' ^^1, 2^ 
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^1,2 2 




^vn? (r^) / 2 


(2.46) 


and, 


C3 I! - 
1,2^ " 


G 1 G 

1,2 p 


1,2 


^ 1,2 ^^ 1,2 


(2.83) 




1^1, 2* ‘^^5,2 ■ ^11,2 ’ ^‘li,2'‘' 3 ^1,2 ^'^1,2 “'^1,2^ 


■ 1,2 


1,2 


R 


- 1 


1,2 


3a 


1.2 


■^ 1 , 2 '"^ 1 , 2 ’^ 1 , 2 


[<^1,2 ^'^1,2^ ““l,2 "^1,2] 


(2.86c) 


Note that the points 1,2 and 8 are not necessarily as shown in Fig. 6-2a. 
Their locations can vary along the J-line depending on the mesh size and the 
Mach number in the vicinity of point 3. Figure 6 -2b presents one example 
when the point m + 2 is an upper boundary point and the point 1’ is outside the 
flow field under consideration. Instead of the point 1^ the point 1 is used in 
the calculation. However, the properties at point 1 are not available nor can 
they be mterpolated because point n + 2 is not yet known. For this particular 
case, the flow properties at point 1 are assumed to be identical to those of 
point m 4* 2 so that the properties at point 3 can be approximated. The same 
method is used to calculate point n + 1. The boundary point n + 2 of the K- 
line can then be calculated without problem. Once the boundary point is deter- 
mined the calculation goes back to point 1 and the right running characteristic 
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m+2 



Upper Boundary- 


Right- Running 
Characteristic Line 


in-2 


Gas Streamlines 



Left -Running 
Characteristic Line 


Normals 


Fig. 6 -2b - Mesh Construction for an Interior Po'int 


inter sectioh-ls made using the two boundary points. Point 1 is then solved 
for and the solution proceeds to the plume boundary. 


ship: 


Step 4: The particle density at point 3K is calculated using the relation- 

j _ 1 

^ U^ (y^^ y^j^^ - y I y^* ) - v| + Ts* ) (^1 " ^3 ^ + ^ 3 ^ 

p1 H -^1 -4) - 

+ 4^44^ “ 44^^ "4 [(4^+4^^^4 " 4^ ■ ^4 4 ^^4~ 4^] 
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Step 5; The particle streamline intersections with the J data line are 
now located (point 8 of Fig. 6- 2a) using the particle streamline relationship; 


J ,J 


4 - ^8 


Step 6; The gas and particle properties at point 8J are now found by 
employing an interpolation scheme as described in Step 3. 

Step 7; The coefficients for the particle compatibility relations are 
now calculated from Eq. (3 -150a). 


^ 8,3 " ^8 


3 ) 
hi' aJ 


C3J , = 


^ 8,3 -^ 8 , 3 / 


- 3 H,3 ^8,3 <H,3 ■ '^8,3^ + ,„j 

^ 8,3 "^ 8,3 


*^8,3'-^8,3' “8,3 ^8,3 


Step 8; The new particle properties at point 3K are then found by ex- 
panding the difference Eqs. (3.148), (3.149) and (3.150) to obtain 


“S “8 ^ ^8, 3 ^^8-3 
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H ~ ^ ^8,3 ^^8-3 


H " ^ 8 / ^.'^ 8 , 3 '^^-3 


with the remaining properties defined as 


and, 


Tj = f(h)^ tabulated 


j 


= f(Re^) tabulated 


0^3 = f(Re|, Pr^) Ref. 12 


9 3*^3 

"^3 ^^ 3 ) 


C 

3 P 3 




(2.46) 


(2.83) 


Q . 1 U 3 * ^^3^ -q3^*^q3^ + 3 C^(T^-T 2 ) 

P 3 ' 3 I 


.j j ' j 

Ai mi ri I 


‘3* 3 3 


(2.86c) 
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step 9; The coefficients for the gas streamline compatibility relations 
are now calculated 





“^5,3 



AL 


5-3 


^ 5,3 "^ 5,3 * 15 , 3 


and 


NP 


5,3 1=1 


/^5/3 ^5,3 [^'^5,3 ■ '^5,3^ '^°®®5,3 ^^^5,3 " '^5,3^ , 3 j ^5 - 3 


,1 


Step 10; The rate of change of entropy in the direction of the gas stream- 
line and the entropy at the point 3K are calculated in the following manner; 


^^5-3 " % 


(3.143) 


and the entropy at point 3K is 


S3 = Sj + Cj 


Step 1 1 ; The rate of change of enthalpy in the direction of the gas stream- 
line and the enthalpy at the point 3K are calculated in the following manner: 


'^^ 5-3 “ ’^ 5,3 '^^ 5-3 " ^3 

' r 


(3.142) 
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and the enthalpy at point 3K is 


Hj = Hg + Tg 3 AS5 2 - Cj . 

P 

Step 12; The rexnaining gas properties to be calculated at point 3K are 
the velocity and flow angle. Expanding equation (3.146) in terms of both the 
right and left-running characteristics and solving for the velocity gives 


01-02 + Qiqi + - B 1 -B 2 + CIi(H3-Hi) + Gl2(H3-H2) + Gi + G 2 - Clj - CI 2 

The new flow angle is then 


Step 13; - The properties T^, p^, u^, , Pr 3 > A^, 


and B 3 are recal- 


culated to reflect the new gas properties. The solution is then checked to 
determine if the desired accuracy has been obtained. If the solution has not 
converged, return to step 2 and repeat the process. The flow of the numerical 
solution is summarized in Chart 6-1. 
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Initial Guess for 
Properties at Interior Point 


Set Up Logic for 
New Point 


Determine Physical 
Location of Interior Point 


Output 


Determine Left-Running and 
Right -Running Characteristic 
Intersection with J-Line 


Check for Desired 
Convergence 


Determine Gas and Particle 
Properties at Intersection 


Calculate New Properties 
at Interior Point 


Determine Particle Density 
at Interior Point 


Calculate Velocity and Flow 
Angle at Interior Point 
Using Compatibility Equations 


Determine Particle Stream- 
line Intersection with J-Line 


Determine AS and AH 
Along Gas Streamline 


Determine Gas and Particle 
P roperties at Inters ection 


Determine Particle 
Properties at Interior 
Point Using the Compatibility 
Equations 


Chart 6-1 - Flow Chart of the Calculation Procedure for an Interior Point Solution 
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6.2 LOWER WALL POINT 


Figure 6-3 presents the mesh construction required to solve for gas and 
particle properties at any lower wall point in Figure 6-1, 


Normal 


Particle 

Streamlines 



Right- Running 
Characteristic Line 


Gas Streamline 
Along Lower Wall 


Fig. 6-3 - Mesh Construction for a Lower Wall Point 


The calculational procedure followed to determine the flow properties 
at a lower all point is very similar to that of an interior point solution. The 
lower wall point is always the first point cadculated on a K-line, Before the 
point 3K can be located, jpoint IJ must first be selected as a function of the 
desired mesh size. Once point IJ is located and its properties determined 
using the interpolation scheme employed for an interior point; the physical 
location (x^, y^) of poinbSK is found by intersecting the right- running charac- 
teristic from point I J and the projection of the streamline from the point 5J, 
The particle density at point 3 K is then calculated using the relationship 


P3 = 


P^5 - ^P\ '^5^ 


(3.163) 


ui 


and is based on the assumptions that: (1) the lower boundary runs parallel to 
the X-axis; and (2) the particle streamlines run parallel to the lower boundary 
when in the vicinity of the lower boundary. 
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Finally, since the flow angle along the lower wall is known, only the right 
running characteristic line is needed to calculate the velocity at the lower wall 
point. Solving for velocity in Eq. (3.146) gives 


[ej-Gj+Q^qj -Bj + CIj(H3-Hj) + Gj - ClJ 



Once the velocity is found, the calculational procedure is continued as 
in the case of an interior point. 

6.3 UPPER WALL, POINT 

Figure 6-4 presents the mesh construction required to solve for the gas" 
and particle properties at any upper wall point in Fig. 6-1. 



The calculational procedure followed to determine the flow properties at 
an upper wall point is almost identical to that of an interior point solution. The 
exception is that the, flow angle on the upper wall is known and therefore, only 

the left-running characteristic line is needed to calculate the velocity at the 
upper boundary point. Solving for velocity in Eq. (3.146) gives 
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*13 = 


03-02 + Q2q2 - B2 + Cl2(H3-H2) + G2-Cl2 


Q, 


6.4 FREE BOUNDARY POINT 


Figure 6-5 presents the mesh construction required to solve for the gas 
and particle properties at a free boundary point. Note that the mesh construc- 
tion is the same as that for an upper wall point. 



Fig. 6-5 - Mesh Construction for a Free Boundary Point 


The calculational procedure followed to determine the flow properties at 
a free boundary point is like that of an upper wall point solution. The one dif- 
ference is that the velocity and not the flow angle at the free boundary point is 
known. Applying Eq. (3.146) to the left- running characteristic and solving for 
the flow angle at the free boundary point gives 


The velocity at the free boundary, though assumed to be known* is not 
easily obtained. The velocity at the free boundary must be solved for itera- 
tively using the relation 


'3 = 


2yRT 

y-1 


(f) ^ - 


1/2 
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An interative solution is required because the gas properties R, T and P 
are functions of a variable 7 . K an ideal gas was being considered, y would be 
constant and v^ may be solved for directly. . 

6.5 PARTICLE LIMITING STREAMLINE PC&SfT 

I" 

Figure 6-6 presents the mesh construction required to solve for the gas 
and particle properties at any particle limiting streamline point in Figure 6-1. 



Fig. 6-6 - Mesh Construction for a Particle Limiting Streamline Point 

The calculational procedure followed to determine the flow properties at 
a particle limiting streamline point is somewhat different from that of an inter- 
ior point solution. The differ ences in the calculational procedure are brought 
about by the fact that the J and K normals are constructed between a gas stream 
line and a particle limiting streamline rather than between two gas stream- 
lines . 


The physical location (x^, y^) of point 3K is found by intersecting the 
normal from the point 7K with the projection of the particle limiting stream - 
line of species J from the point 5J. Once point 3K is located, the gas streamline 
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inters ecUon with the J data line must be located (point 9 of Figure 6-6) using 
the gas streamline relationship 



( 3 , 41 ) 


The properties at this point are found by interpolation between the known points 
4J and 5J. 


The particle density at point 3K for all particle species except particle 
species J is determined in the same manner as for an interior point. The par- 
ticle density of species J is calculated using the relationship 


3,5 6 , , 6 . 6 ,, , 

- ^ - ^ 2 ^ 2 ^ - ^3 ^5^3 + y2) (^3 -x'2> 

(kg - k^) 

and is based on the fact that the entire mass of particle species J is contained 
in the stream tube formed by the particle limiting streamline for species J. 

Once the particle densities have been found, the calculation procedure is 
coptinued as in the case of an interior point. 


+ P: 


u. 


<^54 - y3J^3> - ^2 (4 +>2^ (^5 v ^2^ ^^3 " ^ 2 )] 


(3 


Pc 


h 


^5^ 


5 . / 5 , 6 ^ 
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6.6 PARTICLE LIMITING STREAMLINE-BOUNDARY INTERSECTION 

Figure 6-7 presents the mesh Gipnstruction required to solve for the gas 
and particle properties at the particle limiting streamline -boundary intersection. 



Assume that the particle limiting streamline point (n, K) has just been 
completed and the next point to be calculated is that of a free boundary point 
(n+ 1, K). Under these circumstances, the free boundary point is calculated 
in the normal manner, and then upon completion of the free boundary point a 
comparison is made of the physical location of the two points. If y ^ y » 

an interseGtion of the particle limiting streamline and free boundary is indi- 
cated. The physical location of {n + 2, K) is found by intersecting the projec- 
tion Pf the streamline from the point (m+ 1, J) and the particle limiting stream- 
line between the points (m, J) and (n, K). The gas and particle properties at the 
point (n+ 2, K) are then found by interpolation between the points (m, J) and (n, K). 

Once the properties at the point (n+ 2, K) are found, the point (n+ 1, K) is 
assumed to no longer exist. The point (n+ 2, K) is used as the base point the 
next time an upper boundary point is calculated. 

The description of the entire flow field now proceeds onto the calculation 
of the lower boundary point of the new line. 
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6.7 SHOCK WAVES 


Unlike the interior point solution, there are two points assigned to every 
shock wave point. To determine if the shock wave is either left-running or 
right- running, one can check the shock points (m, J) and (m-1, J) in Figure 6-8. 

If the point (m, J) is a shock upstream point and the point (m-1, J) is a shock 
downstream point, the shock wave is a left running one as shown in Figure 6- 8(a). 
If the case is reversed, the shock wave is right- running as shown in Figure 
6-8(b). Once the shock wave is detected, the new shock point (n, K) is found by 
extending the shock wave from point (m, J) and the normal from point (n-1, K). 
Generally, first, the properties of the shock upstream point are found according 
to the upstream reference properties and then the properties of the shock down- 
stream point are calculated with the oblique shock relations using Eq. (4.22). 

The velocity of the shock downstream point is also calculated by using the 
shock downstream reference properties with one characteristic. The velocity 
calculated by both methods is compared and the shock angle is adjusted until 



Right- Running 


Characteristics 


50 n-l 


Left -Running 
^Characteristics 



Right- R unning 
Characteristics 


K 

Normals 


Normals 


Fig.6-8a - Mesh Construction for a 
Left-Running Shock Wave 
Point 


i 


Left-Running 

Characteristics 


Fig . 6 - 8b - Mesh Construction for a 

Right-Running Shock Wave 
Point 
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the velocity calculated by both methods satisfies the convergence criteria, 

(Note that the streamline is generally terminated in front of the shock wave 
and a new one started from the shock downstream point.) 

6.7/1 Left-Running Shock 

As shown in Figure 6-8(a), the shock point location is found first and the 

I shock upstream point (n+ 1, K) is treated as an interior point. The calculation 

i ■ '■ ■■ \ 

is carried out as described in the preceding paragraph. 

I ,' ■ ■ ■ ■ ' 

* Under normal conditions, the calculation for point (n-1, K) in Figure 

6- 8(a) is slightly different from that for the regular interior point. As can be 
seen; the right-running characteristic from that point intersects the shock wave 
before intersecting the J-line. The base point properties on the right- running 
characteristic are interpolated between point (m-1, J) and (n, K). Because of 
this, points (n-1, K), (n, K) and (n+ 1, K) are handled simultaneously. Under 

I certain conditions, the right- running characteristic from point {n-2, K) must 

i H . . 

;\be added to this group if its right -running characteristic intersects the shock 

I wave before intersecting the J-line. 

I ■ ■ 

6.7. 2 Right- Running Shock 

1 ' 

* . Assume that a right-running shock is detected when the point (n ,K) as 

? shown in Fig . 6-8(b)^^-i's ECtFempted. The calculation for this point is termi.- 

* nated and the intersection of the shock wave with the K -line is found. The 

^ shock upstream point (n,K) is treated as an interior point and the shock points 

I calculation is carried out as usual. Under certain cases, the properties of the 

shock upstream point are interpolated from the two known points on the K -line 
i if the intersection of the shock wave with the K-line falls between these known 

■ V'-"- 

points. For instance, if the point (n , K) in Fig. 6-8(b) is complete, the prop- 
I erties of the shock upstream point (n, K) are interpolated from points (n - 1 , K) 

^ and (n^, K) (see Section 6.8.1 Upper Compression Corner) . 

i The calculation of the first few interior points on the K-line behind the 

^ right -running shock wave is slightly different from the regular interior point,, 

i ■ 

i ' . . ^ ' 
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As shown in Figure 6-8(c), the left- running characteristic from point (n+2, K) 
intersects the shock wave and its base point properties are interpolated be- 
tween points (rn, J) and (n+1, K). The left-running characteristic from point 
(n+3, K) would have intersected the shock wave but the base point properties 
are interpolated between points (m, J) and (n+2, K). 


Normals 


Right- Running 



Fig. 6 -8c - Mesh Construction for an Interior Point 
Behind a Right -Running Shock Wave 


6.8 COMPRESSION CORNER POINT 

The locations of compression corners in the flow field are defined by the 
boundary conditions of the problem. Whenever a boundary point is calculated, 
the location of this point is always compared to the limit of the current boundary 
equation. If the boundary point is already beyond the limit, the properties of the 
limiting point are interpolated between the boundary point on the J-line and the 
newly calculated boundary point whose location has exceeded the limit of the 
current boundary. The calculation thereafter will depend on the type of dis- 
continuity and the type of following boundary. 
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6.8,1 Upper Compression Corner 


Figure 6 -9a shows a typical compression corner at the upper boundary. 
The interior point K) and the upper boundary point (n4l, K) are first calcu- 



Fig.6“9a - Mesh Construction for an Upper Fig.6-9b - Mesh Construction for 

Compression Corner Point a Lower Compression 

Corner Point 

lated without any knowledge of the existence of the compression discontinuity. 

Then the properties at the limiting point (m+3, J) are found by interpolation be- 
tween the points (m+2, J) and (ri^>l,K). Since the boundaries are already given, 
the* amount of turn the flow must make can be found directly from these bound- 
aries. With the information of the known turning angle and the flow properties 
at the corner, point (m+3, J), a right-running shock is defined and the properties 
of the shock downstream point (m+4, J) can be calculated with the oblique shock 
relations using Eq. (4,22) . After the shock wave is established^ it is extended 
from the discontinuity corner to intersect the just calculated K -line. The flow 
properties at the interseGtion point (n,K) in Fig. 6-9ci, are interpolated between 
points (n-l,K) and (n'jK) which are air eady calculated. The properties 
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thus obtained are for the shock upstream point (n, K). The calculation for the 
shock downstream point and the boundary point is as follows. The properties 
of the shock downstream point (n+1, K) are first calculated using the oblique 
shock relations. A normal is then drawn from point (n+l,K) to Intersect the new 
boundary at point (n+2, K). The flow properties at point {n+2, K) are at first 
assumed to be the same as those of point (m+4, J). Note that, by doing so, the 
entropy of point (n+2, K) is automatically set equal to that at point (m+4, J). 

The velocity at point (n+2, K) is then calculated by using the left- running charac- 
teristic with its base point properties interpolated between points (m+4, J) and 
(n+1, K). The other properties at point (n+2, K) are also adjusted accordingly. 

The velocity at point (n+1, K) is then calculated by using the right- running 
characteristic with its base point properties interpolated between points (m+4, J) 
and (n+2, K). The flow velocity thus calculated is compared with that calculated 
from the oblique shock relations based on the shock upstream properties. The 
shock angle at the K-line is adjusted and the calculation repeated until the velocity 
calculated with both methods satisfies the convergence criteria. At the begin- 
ning of every new calculation cycle, point (n, K) is relocated with the new av- 
erage shock angle between the points (m+3, J) and (n, K), and the new average 
normal between the points (n-1, K) and (n, K). The average normal is also 
used to locate point (n+2, K). 

6.8.2 Lower Compression Corner 

Figure 6-9b shows a typical compression corner at the lower boundary. 
Assume that the J' -line has just been calculated. The lower boundary point 
(1' , K' ) is first calculated without any knowledge of the existence of the com- 
pression discontinuity. The flow properties at the compression corner are then 
calculated in the same manner as for the upper compre ssion corner and a left- 
running shock wave is established. A new J-line in the shock upstream region 
is created emanating from the lower compression corner. 

To start a new K-line in the shock downstream region corresponding to the 
new J-line in the shock upstream region , the following method is adopted. The 
intersection of the shockwave from point (1, J) and the streamline (if the shock 
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angle is very small, the right- running characteristic is used) from point (2, J) 
is found first. This point, the shock upstream point (3, K), is treated as an 
interior point and the oblique shock relations are used to calculate the pro- 
perties of the shock downstream point (2, K). An average normal is drawn 
from this point to intersect the new lower boundary, which is designated as 
(1, K) in Figure 6-9(b). The properties at point (1, K) are initially set equal 
to those at point (m+4, J); its velocity is recalculated by using the right- run- 
ning characteristic with its base point properties interpolated between points 
(2, K) and (m+4, J). The other properties at point (1, K) are adjusted accord- 
ingly. The velocity at point (2, K) is then calculated with the aid of the left- 
running characteristic with its base point properties interpolated between 
points (m+4, J) and (1, K). The rest of the calculation is carried out in the 
same way as in the case for the upper compression corner. 

6.9 EXPANSION CORNER - PRANDTL-MEYER FAN 

Expansion corners are found when the slope of the solid boundary has a 
discontinuity in such a way that the flow must negotiate a sharp turn away from 
the oncoming flow or when the flow issuing into an open space is under-expan- 
ded in a nozzle or a flow channel. The expansion takes place eithei* in a chan- 
nel or at the channel exit. The amount of turn the flow must make and the fi- 
nal flow velocity can be found with the aid of the Prandtl- Meyer expansion re- 
lation, Eq. 5-1, and the boundary conditions. 

After the total turning angle is found, it is subdivided into a number of 
small turning angles. The flow properties at the expansion comer correspon- 
ding to each of the small turns are calculated and the angle' s corresponding 
characteristic lines emanating from the corner are constructed. The final 
line of the expansion is a streamline rather than a characteristic line and, 
normally, the Mach angle corresponding to this streamline is substantially 
greater than the small turning angle. An extra point is inserted to fill the 
gap after the expansion corner computation has been completed. Note that 
this final expansion line, a streamline, can be a solid boundary or a free > 
boundary. 
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In the example shown in Figure 6-10, the turning angle, is subdivided into 
four small turns. The corresponding characteristic lines emanating from the 
corner are designated as A, B, C and D in Figure 6-10. The final expansion 
line E is a streamline corresponding to the last corner point (d, J). The pro- 
perties at the corner corresponding to every small turn are stored in (a, J), 

(b, J), (c, J) and (d, J) arrays. The properties at the corner before the ex- 
pansion takes place can be found as usual and are stored in (m, J) array; its 
characteristic line is designated as M in Figure 6-10. 


A weak shock is initialized point (n+4, K) for the upper corner expan- 
s\o!h point or (n-4, K) for the lower ..corner expansion, which is done simply by 
defining a double point at those points with identical flow properties. The initial 
shock angle is assumed equal to the Mach angle. 


6.9.1 Upper Corner Expansion 

As shown in Figure 6- 10 (a), the point (n-1, K) is considered to be com- 
plete. A left-running characteristic is drawn from the point (n- l,K) to inters’ect 
the right-running characteristic from point (m, J) at M, which is an interior 
point. The properties can be calculated in the usual manner. A normal is drawn 
from point (n-1, K) to intersect the right- running characteristic from point (m, J) 
at n. The properties at this point are then interpolated betw'een points (m, J) and 
M. Similarly, the properties at points (n+l, K) through (n+4, K) can be calcu- 
lated. Point (n+5i K) is calculated by the same technique provided that the point 
E is treated as an upper boundary point. If the distance between points (n+4, K) 
and (n+5, K) is greater than other meshes in the expansion fan. One extra point 
(n+6, K) is inserted and its properties are found by interpolation between points 
(n+4, K) and (n+5, K). jf 


6.9. L Lower Corner Expansion 

This case is shown in Figure 6- 10(b). When an expansion at the lower 
boundary is detected, the normal line starting from the corner point (m, J) is 
completed first. The flow properties at the corner and the corresponding 
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characteristic line for each small turn can be found by the same method de- 
scribed in the preceding paragraph. Point (n, K) is the intersection of the 
right- running characteristic from point (m+1, J) and the left- winning charac- 
teristic from point (m, J); its properties can be calculated without problem. 

The flow properties at point A can also be calculated. A normal from point 
(n, K) and the left- running characteristic from point (a, J) intersects at (n-1, K); 
flow properties at this point are then interpolated between points (a, J) and A. 
Points (n-2, K) through (n-4, K) can be calculated by the same method. Point 
(n-5, K) is a lower boundary point which can be calculated with the technique 
used to find the point (n+5, K) in Figure 6-10(a). An extra point may also be 
inserted between points (n-4, K) and (n-5, K) if needed. 

6.10 SliiP LINE POINT 

Sliplines are normally found behind shock wave intersections and in mul- 
tiple streamtube flow where each streamtube retains its own identity. For every 
slipline in the flow field, two points are assigned and initialized at the input line. 
The conditions which the slipline points have to satisfy are that the flow angle 
and the pressure on both sides of the slipline must be the same. 

The point (n, K) in Figure 6- 11 (a) is found by extending the slipline from 
point (m, J) and the normal from point (n-T, K) as usual. First, the flow angle 
at the new slipline points is assumed to be identical to that of point (rn, J), and 
the velocity at points (n, K) and (n+1, K) are evaluated by using II- and I-charac- 
teristic respectively. Then, the pressure at both sides of the slipline is calcu- 
lated and compared. If a pressure difference exists, the flow angle is adjusted 
back and forth until the pressure is balanced on both sides of the slipline. 



Fig, 6- 11a - Mesh Construction for 
a Slip Line Point 

■■ rf- 


Fig.6-llb - Mesh Construction for 
a Slip Line Point Where 
a Characteristic Inter- 
sects a Slip Line 


During the calculation, crossings of the characteristic line with the slip- 
line are not allowed. For example in Figure 6- 11(b), the I- characteristic of 
point (n, K) intersects the slipline first before intersecting the J-line. The 
reference properties for the calculation of point (n, K) on I- characteristic are 
interpolated between points (m+1, J) and (n+1, K) if point (n+1, K) is known; when 
the point (n+1, K) is still unknown, the properties at point (m+1, J) are used. 

The slipline points (n+1, K), (n+2, K) and (n+3, K), (n+4, K) can be handled as 
described in the preceding paragraph. For the case shown in Figure 6- 11(b), 
the points (n, K) (n+1, K), (n+2, K), (n+3, K) and (n+4, K) are iterated as a group, 
since the calculation of any point in this group is interrelated. 
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Section 7 

NUMERICAL INTEGRATION OF THE CONSERVATION EQUATIONS 

Examination of the mass flow rate, momentum flux and energy flux 
through a nozzle and exhaust plume will yield considerable information about 
the nozzle performance and the status of the flowfield numerical solution. 
The momentum flux is used to determine the thrust produced by the nozzle 
while the thrust ratioed to the mass flow rate is a measure of the nozzle 
performance. Deviations of the fluxes in the conservative parameters from 
surface to surface (Fig. 7-1) are indicators of the validity of a numerical 
solution such as the one used in the gas -particle code. 



Fig. 7-1 - Schematic of Conservation Parameters in a Nozzle Plume 
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7.1 CONTINUITY 


The mass flow rate through a surface is calculated by numerically inte- 
grating the incremental mass flowthrough adjacent streamtubes. A streamtube 
is defined as the flow area bounded by adjacent data points (Fig.,;7-2a) and the 
incremental flow rate is calculated as follows. . 


r 



Fig. 7-2c - Schematic of Surface 
Normal Unit Vector 

Fig. 7-2 - Schematic of Elemental Geometry on a Data Surface 


Let points a and b be any two points on a data surface within the flow field. 
The incremental flow area for axisymmetric or two-dimensional flow is given 

by ■ 

dA - (Ztt r) dL . (7.1) 




■ I 
iJ 









Wst 


i 
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In terms of the coordinate directions, the incremental area is 

A ? 7 >1/2 

dA = i2nr)° {(dr)'^ + (dx)^| 

The flow area in the direction normal to the surface, dL, is now 

dA = (2ir r)^ ] (dr)^ + (dx)^l 


(7.2) 


(7.3) 


where the unit normal vector, ©^, is defined by the geometry of Fig. 7-2c as 

A 


®n ~ cose/? i + sin<^ j ; 


(7.4) 


and the angular location of the unit normal vector relative to the x-axis is 


(p = - tan 


'») 


(7.5) 


The quantity of mass flowing through the incremental area is found from 
the relation 

dm = p|^|dAeq*e^; (7.6) 

where the unit velocity vector, 0 is defined by the geometry of Fig, 7-2b. 

^ .. 

The relation for m is now 


• I I t A. I ^ ^ \ 

m = p I q |dA |cos9 i + sin9. j !• • (cos(/> i + sirup j j 

Performing the indicated dot product the relation reduces to 

^ = P [‘I i d-A’ |cos9 cos(^ + sin9 sin<p| . 

The trig ometric identity for sine -cosine function is 

cos9 cos(p + sinS sirup = cos . 
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Substituting, the incremental mass flow is 


dm - pj q [dA co8{Q-cp) 


(7.7) 


The total mass flow through the surface is simply 


m = dm 


(7.8) 


for all the incremental flow areas. 

So far the discussion has been directed toward a general expression 
for the mass flow rate. For gas -particle flows the surface geometry is the 
same. The only change is in the definition of the flow variables. The ex- 
pression for the respective phases are as follows: 

Gas Phase 


m = p q dA cos(0 -(p) ; 


(7.9a) 


Particle particle) 


= P”' 1^ dA cos (9^ -<^P) 


(7.9b) 


Mixtur e 


7.2 MOMENTUM 


NP 

mm = ^ 


J 


j=l 


(7.9c) 


As with the mass flow relations, a general expression for the momentum 
flux will be developed. This expression will then be extended to include the 
treatment of gas -particle flows. 
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REPRODUCIBILirY OF THE 
OMINAL PAGE IS POOR 


The general expression for the momentum flux is 


Mm = jP6+pqqj*dA 


(7.10) 


where 6 is the Kroneker delta and dA is given by Eq, (7.3). Substituting for 
dA and performing the indicated algebraic operations, Eq. (7.10) becomes 


+ P q q . n 


dA 


(7.11) 


where dA is now defined by Eq. (7.2). Equation (7.11) yields an expression 
for each coordinate direction in the equation in terms of the scalar compo- 
nents given by the expressions 


Mm = |p{cos(/?i + sin(^j) + p(ui + vi) (ui* + vj) • (cos</^ + sin<^j)|. dA ; 


M 


M 






|^P(cos</?i. + sinCpj) + p(u ij i' uvlj + uyji + v ji) • {cosipi + sin(j^?j) | dA ; 


^ A 2 A A A 2 A ) 

Mm - yP(cosc/?i + sin^j) + p(u cos^ + uv coscpj + uvsinr^i + v sin^j)} dA 


The scalar equations in the coordinate directions are then 

Mm = coscp + p (u^ coscp + uv sin(p)| dA , 


and 


M 


( 2 ) 

M = |P sin</> + p(uv coscp + v sin(/?) dA; 
yy i 


Substituting for the velocity components 


Mm = [P coscp +pq^ (cos9 cos9 cos(iC? + cos9 sin9 sin^)| dA 


and 


M, , = I p smU) +Pq^ (cos9 sinO coscp + sinO ^'inO sirup)] dA 

M_ ■ ( ) 
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n 

U 


n 


Factoring common terms and substituting the trigometric identity as 
the final result for the incremental momentum'-flvix is 

before 

y 

< 

II 

2 

P cos<p + pq cos(9 - <^) COS0 

dA , 

(7.12a) 

li 

and 




J 

= 

Mr 

2 

■Psin^ + pq cos (Q -if) sin9 

dA 

(7.12b) 

WsV* 

ll 


The only difference between the expression for the gas and particle 
contributions to the momentum flux is that the particles do not contribute to 
the pressure of the mixture; The respective expressions for the momentum 
flux through the surface are; 


A 


Gas 

= |p cos^ + pq^ cos (9 -(/?) cos9 

X 

= |p sin(^ + pq^ cos(9 -(/?) sin9| A 

Particles (y^ Particle) 

= p^ (q^) cos{Q -(p) cos(9'^') A 

X 

• * O * 

= p^ (q ) cos(9 -(/?) sin(0'^) A 


(7.13a) 

(7.13b) 


(7.13c) 

(7.13d) 


It 


fi i 

(1 IS 






h 

f p 


ri 


Mixture 


M 


M 




m 

X 


m 


NP 

= M + V' 

M Z-/ 

X X 


= M, 


NP . 

j = I 


(7.14a) 


(7.14b) 


rii 

n c 


/r * . 

k 




7“6 


LOCKHEED - HUNTSVILLE RESEARCH A ENGINEERING CENTER 


7.3 ENERGY 


The general relation is 

E = (h +|-q^) m (7.15) 

Applied to gas and particulate phases the respective relations for a given 
surface are: 

Gas ’ 

E = (h + -I- q^) m 
th. 

Particles (j Particle) 

E^ = (hj + - (q^)^) 

Mixture 

NP , 

E^ = E + 2]] (7.16c) 

j=l 

7.4 NOZZLE PERFORMANCE 

Performance parameters considered are the thrust produced and the 
associated specific impulse, Igp* The force vector is calculated as the 
mixture momentum on the initial data surface (Fig. 7-1) plus the pressure 
force on the nozzle wall. Mathematically this is expressed as 


(7.16a) 


(7.16b) 


= % + /Pw <■'•1^) 

The integral term in Eq. (7.17) is obtained by numerically integrating the 
pressure force along the nozzle wall. Consider the elemental surface, dL, 
bounded by the nozzle wall points a and b on any two successive data 
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Data Point, b 


Nozzle Wall 


Data Point, a 



Note: is the Nozzle 

Thrust 

Fj.=0 for axisym- 
metric Flow 


Fig. 7--3 - Schematic of Pressure Force Acting on the Nozzle Wall 


surfaces. The slope of the surface, dL, is given by 


(f) • 


( 7 . 18 ) 


However, the force due to the pressure loading acts normal to the surface 
which is given by the angular relation 90 + cp. The local surface area over 
which the pressure is acting is 


dA = {2w rf ! (dr)^ + (dx) ^ l 


I' 


I 


( 7 . 19 ) 


and the unit normal to the surface is 




(> A 

sin(^,L + cosi;'7j , 


(7.20) 
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so that 


dA = dA n 

1/2 


dA = + (dx)^ 


A A 

sin<^i + cos((p)j 


(7.21) 


Equation (7.21) can be generalized to reflect the loading increment to either 
an upper or lower wall in the following manner. Define a parameter A which 
has a value of + 1 for an upper or lower wall, respectively. The orientation 
of the surface unit normal with respect to the x-axis is 


^1 "" ^ 


From the sine and cosine trigometric identities Eq. (7.22) reduces to 


and 


cosXPy - - sincp 


sin(p^ = coscp 


so that the unit normal vector is now given by 

A A 

n = coscpj^i + sin^j^j ; 

, li ■ 

and the surface area term is now If 


dA 


w 


= |(dr)^ + (dx)^l |cos((/?j)i + sin(^j)j| 


(7.22) 


The resultant force relation for the nozzle is then 

NS 1/2 

^ + 5^ jp (2xr)^ (dr)^ + (dx)^| {cos(<??i)i + suK^Jj] (7.23) 

^1=1 V 

and the nozzle thrust is just F^. The specific impulse, is then 
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(7.24) 


I = F /m_ (7.24) 

sp x' rrij. 

7.5 PERCENT CHANGES IN THE CONSERVATION QUANTITIES 

One check on the validity of the numerical solution is to monitor the 
change in the conservation quantities mass flow, momentum and energy. If 
these quantities do not change relative to the initial surface values, then the 
mass, momentum and energy have been conserved; this is a good indication 
that the numerical solution is proceeding satisfactorily. However, since a 
numerical solution is employed some changes in the conservation quantities 
is to be expected. The objective then is to keep the error as small as possible. 
When large errors occur then the input data, wall equations, mesh construction, 
etc., should be checked. 

The changes in the conservation quantities are found by comparing the 
integrated value at a given data surface relative to the initial data surface. 
These are as follows: 

Mass Flow Rate (Mixture): 


%Am =: 


(rrij^ - riij) . 100 


(7.25) 


Momentv^m (Mixture) 


%AM = 


(|^| - |Mj )|. 100 


(7.26) 


%AE = 


{E^ - Ej) • 100 


(7.27) 


The individual changes for the respective phases are computed in a similar 
fashion as specified by Eqs.’ (7,25) through (7.27). One important item to 
remembef''is';;.that when gas particle flows are being considered the 
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solution is coupled in the sense that momentum and energy is exchanged 
between the phases. This means that relative large changes in the momentum 
and energy quantities for each phase can occur. What should be of concern 
for gas -particle flows is the change in the mixture conservation quantities 
(momentum and energy). The respective mass flows should be conserved 
since the assumption has been made that the phases do not exchange mass 
via chemical reactions, condensation, etc. 

Finally, the specific impulse for a nozzle can be calculated in one of 
several ways. One way has already been outlined where the thrust is computed 
as the sum of the initial surface momentum and the pressure force acting on 
the nozzle wall. The specific impulse is then calculated by dividing the thrust 
by the mass flow rate through the initial surface. Using the initial surface 
mass flow rate is a matter of choice and can be argued, from either a pro or 
con view. However, use of the initial surface mass flow rate eliminates the 
use of one variable which would be numerically calculated. Computation of 
the specific impulse by the second procedure is as follows. 


The thrust generated by the nozzle is first obtained by calculating the 


momentum in the nozzle axial direction (M 


M 


, Eq. (7.14a)) . The specific 


m 


X 


impulse is then calculated by dividing by the mixture mass flow rate. This 

gives a means to compare the I calculation previously discussed. The 

®P 

change in I is 
sp 


(M 


%AI 


, , /m - I ) * 100 
M ' m sp' 
m 

X 


sp 


I 


sp 


(7.28) 


where the term 1^^ is obtained from Eq. (7.24). 


7-11 

LOCKHEED • HUNTSVILLE RESEARCH & ENGINEERING CENTER 


Section 8 

CONCLUDING REMARKS 


The preceding sections have presented a detailed development and 
description of a numerical solution for the treatment of supersonic gas- 
particle flows* A computer code has been developed and is described in 
Volume II of this report. Several options have been included that necessitated 
a somewhat generalized set of input data. These options include; 

• Gas-Particle Flows 

• Chemical 

• Equilibrium 

• Kinetics 

• Chemically Frozen 

• Constant Thermodynamics 

• Single Phase Solution 

• Non-Isoenergetic Flow 

• Non - Continuum Flow 

• Performance Calculations 

A primary consequence of this work is the extension of gas -particle 
solutions to treat chemical kinetics for nozzle-exhaust plume flow fields* 

Since the code has the option of treating single phase flow^ chemical kinetics 
can also be included in liquid propellant motor analyses as well as solid 
propellant calculations. The method by which the kinetic equations are 
modeled also permits thermal nonequilibrium to be treated. 
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In its present form, the code has the capability of producing data for 
the following applications; 

• Gas/Gas -Particle Impingement (Heat Transfer-Loads) 

• Rocket Nozzle Performance (Thrust, 1^^) 

• IR Signatures (Radiating Species) 

• RF Attenuation (Electron Densities) 

• Plume Radiation (Radiative Heat Transfer Gas/Particles) 

• Vehicle Base Pressure '' 

• Base Heating (Convection-Recirculation) 

Future development may be pursued toward developing the capability of handling 
imbedded subsonic regions (Mach disk) and mixing. 
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Appendix A 

SYMBOLS AND NOTATION 


Symbol 


Description 


f 

Kn 

M 

Nu 


drag coefficient, dimensionless 

incompressible drag coefficient defined by Eq. {A-2) 
Knud sen number, dimensionless 
Mach number, dimensionless 
Nusselt number, dimensionless 


Pr 

q 

r 

R 

Re 


Prandtl number, dimensionless 
velocity, ft/sec 
particle radius, ft 
gas constant, ft^/sec^/°R 
Reynolds number, dimensionless 


T temperature, °R 

Greek 


r 

fJ- 

P 

Subscripts 

C 

FM 

I 


Superscript 

j 


specific heat ratio, dimensionless 
viscosity, Ibf-sec/ft^ 
density, slug/ft^ 

indicates the quantity pertains to the continuum flow regime 
indicates the quantity pertains to the free molecular flow regime 
indicates the quantity pertains to incompressible flow 


indicates the quantity pertains to a particle species 
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APPENDIX A 


Particle drag and heat transfer quantities are calculated in terms 
of non-dimensional coefficients. These coefficients are functions of 
the local Reynolds number or the Knudsen number depending on the 
parameter being calculated. The following discussion summarizes the 
calculations employed in the gas-particle flow analyses. Details of 
the theory are given in the cited references. 

A.l PARTICLE DRAG COEFFICIENTS 

Particles in solid rocket motor nozzle and exhaust plume flow- 
fields encounter a wide range of flow regimes. Consequently any 
expression chosen for the drag coefficient must have the ability to 
take this into account. In a rocket nozzle and exhaust plume flow field 
the particles are continually being accelerated by the turbulent flow which 
has a temperature different than the particles. This situation is further 
complicated by rarefaction effects because the mean free path in the 
gas is comparable to dimensions of the particle boundary layer. 

Obviously the **standard drag curve” of a sphere as a function of 
Reynolds number (Ref. A-1) does not fit this description. This drag 
coefficient can only be applied to 

• Single solid sphere 

• Constant velocity 

• Still, isothermal, incompressible flow of infinite extent 

In early two-phase studies the sphere drag coefficient was taken 
as the Stokes drag coefficient; 
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This expression, however, is valid only for Re L 1 if rarefaction 
effects are not taken into account. 

The Stokes drag coefficient is now used as a reference drag co- 
efficient. An incompressible drag parameter defined by the relation 

fj = Cj,VCDLke = ’ NP <A-2) 


was tabulated as a function of Reynolds number by Kliegel (Ref. A-2), 
This was corrected for rarefaction effects (Ref. A- 3) by the relation 


C D = C 


Dt 


(1+ 7.5 Kn) (1+2 Kn) -1 1.91 Kn^ 

(l + 7.5Kn) (l + 3Kn)+ (2.29+ 5.16 Kn) Kn' 


} 


(A- 3) 


The Knudsen number , Kn, (a measure of flow rarefaction) is related to 
the Mach and Reynolds numbers by 


Knj = 1.26 MVReJ . 


(A-4) 


= 1 Aqj|/vyRT . 


(A- 5) 


A more recent formulation fc the drag coefficient is taken from the 
work of Crowe (Ref. A-4). The drt'.g coefficient is defined as 
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i i 



^Di : 


where 

Kn(C^ - .48) 

T, .7^ inc 

-Kn © g~g Re 

A.2 PARTICLE HEAT TRANSFER COEFFICIENTS 

At present there appears to be no specific experimental heat trans- 
fer data for spheres of the particle sizes in the flow regimes encountered 
in a rocket nozzle- exhaust plume flow field. A review of semi- empirical 
expressions is given by Miller and Barrington (Ref. A-5). Until further 
work is done that is directly applicable to particle flows the continuum 
expression of Drake (Ref. A-6) 

Nu^^ont "" ^ / 

vvhere 


= 


i 


1.1 

+ 1.1 (Kn) 


-.3 


exp 


}{■ 


exp 


Re'^ = 2p I Aq | rV fl 


modified for rarefaction, effects by Kawanau (Ref. A-6) 


Nu 


Nu = 


cont 


1.0+ 3.42 M/Re'^ Pr Nu 

cont 


appears to be the most applicable to two -phase flows. 
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Appendix B 

SYMBOLS AND NOTATION 


Symbol 

h 

H 

P 

q 

s 

T 

Greek 

■n 

P 

A 

Subscripts 

o 

f 

Superscript 


Description 

local enthalpy, ft^/sec^ 
total enthalpy, ft^/sec^ 
pressure, Ibf/ft^ 
velocity, ft/sec 
entropy, ft^/sec^/°R 
temperature, 

weight flow ratio of oxidizer to fuel, dimensionless 
3 

density, slug/ft 

denotes a change in a variable over a step length 

indicates the quantity pertains to the oxidizer 
indicates the quantity pertains to the fuel 

denotes a vector quantity 

denotes an average value over a step length 


FBEGEBIUG r-r A>TTr 
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Appendix B 


The bulk of this document is concerned with the discussion of the 
analysis of gas-particle flows which have been shown to be non-' Lsentropic 
as well as non- Lsoenergetic in the gas phase. However, one of the 
original objectives of the RAMP code development was to provide for a 
single phase (gas) Galculation as an option. Computationally this 
entails simply considering the coefficients reflecting particle contri- 
butions tq the compatibility equation as numerically equal to zero. As 
expected, there is no change in entropy and gas total enthalpy along 
gas streamlines, that is, the flow is seen to be isentropic and iso- 
energetic. However, for many liquid propellant motor applications a 
fuel rich mixture of the chamber combustion products is used to 
film cool the nozzle walls. This results in fuel striations across the 
nozzle where each streamline has a different energy level associated 
with it (i.e. non-isoenergetic flow). 

This case can be treated with a non-equilibrium chemical analysis 
by simply specifying the appropriate initial gas total conditions on each 
streamline. Recall that the pressure-density-velocity form of the 
compatibility equation is then used. Gas chemical concentrations are 
required for each streamline so that the gas thermodynamics are 
computed locally as the solution proceeds downstream in the nozzle- 
exhaust plume flowfield. However, if the flow is to be considered in 
chemical equilibrium or equilibrium/frozen, the applicable compatibility 
relations must x-eflect the non-isoenergetic nature of the striated flow 
condition. These relations are developed in the following paragraphs. 

To begin, the development of the species continuity equation 
(starting with equation 2-6) could be replaced by atomic conservation 
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equations. Moreover if the conservation of those atoms associated 
with the fuel and the oxidizer were considered then, for steady state. 


V * (pft) =0 (B-1) 

and 

V • (Po^) = 0 

would result. 

If the weight flow ratio of oxidizer to fuel (O/f ratio) is denoted by q 
then Eq. (B-l) and B-2) are satisfied if, 

• vn = 0 (B-3) 

and ■■ T. ■ ■ 

V • (pq) =0 (B-4) 

The assumptions inherent in arriving at Eq. (2-50) did not involve 
isoenergetic flow so that the momentum equation remains valid. The • 
energy Eq. (2-105), however, remains valid only along each streanxline 
and must be replaced by 

'-*2 , 

h + I = Hei) (B-5) 

The compatibility equations which apply along each Mach line are unaltered 
by the non-isoenergetic analysis since it is constructed based on the 
momentum equation and the global continuity equation (neither of which 
are altered by the non-isoenergetic flow phenomena). 
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The numerical solution to the governing equations is not greatly- 
affected by the modification to the compatibility equation. In Eq. (3-146) 


AP 


’ 1.2 




1.2 


1,2 \^ 1,2 


^1,2 ^'>1.2 


(B-6) 


while for a lower pressure boundary 


1 ( ^^1 2 — 

<"3-^l,2) <== A— ^ +qi.2^qi,2 

^ 1,2 \ P 1,2 


(B-7) 


The finite difference analog to the streamline equation may be 
written: 

a. interior point 

r,3=q5 

b. an upper wall or upper free boundary point 

ri3=hi 

c. a lower wall or lower pressure boundary point 

Ti3=q2 

The oblique shock relations of Section 4 are unaffected since the 
discussion concerns a single streamline. All references to flow variables 
are understood to be a function of the o/f ratio. Since that ratio does 
not change across the shock the only alteration necessary is that of 
Eq. (4-11) which is altered to read 

^2 “ ^2*^1^^ ^2 ^ P^®2'* ^2' ^1^ (B-8) 
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Since the Prandtl-Meyer calculation of Section 5 is also concerned 
with a single streamline, hence a constant value of it to is unaffected. 

It can readily be seen from the above discussion that the alterations 
necessary to incorporate the non-isoenergetic analysis are straight- 
forward from an analytic point of view. 
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Appendix G 

CHEMICAL EQUILIBRIUM CALCULATIONS 
IN GAS- PARTICLE FLOWS 
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Appendix C 

SYMBOLS AND NOTATION 


Symbols 


A 



c 


P 

H 

M 


w 

N 

N 


P 

P 

q 

R 

s 

T 


u, V 


V 

w 

X 

Greek 

P 

0 


Description 

parameter defined by Eq. (2-46), l/sec 

parameter defined by Eq. (2-103) 

specific heat at constant pressure, ftVsecV°R 
2 2 

total enthalpy, ft /sec 

molecular weight, Ibm/lb mole 

number of moles 

number of discrete particles 
2 

pressure, Ibf/ft 
velocity, ft/sec 
gas constant, ft^/ sec' 
entropy, ft^/ sec‘ V°R 
temperature, 

velocity components, ft/sec 
velocity, ft/sec 
flow rate, Ibm/sec 
position coordinate, ft 

3 

density, slugs/ft 

inclination of the flow vector with respect to the x-axis, deg 
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SYMBOLS AND NOTATION (Cont^d) 


Subscripts 

c 

g.G 

m 

P 

Superscript 

j 


denotes chamber conditions 

indicates the quantity pertains to the gaseous species 

indicates the quantity pertains to the gas-particle 
mixture 

indicates the quantity pertains to the particle species 
!' indicates the quantity pertains to the particle species 
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Appendix C 

One of the basic assumptions made in the development of equations 
describing the gas -particle flow is that there is no mass transfer between 
the phases-, Consequently, calculations treating chemical reactions in gas- 
particle flows must reflect this assumption. Prozan (Ref. C-1) has shown 
for chemical equilibrium that the gas thermochemistry calculations can be 
uncoupled from gas-dynamic considerations. Gas thermodynamic properties 
(Ref. C-2) for chemical equilibrium or frozen flow are precalculated for nozzle 
chamber conditions. Tabulated thermodynamic data are then constructed in 
terms of the independent variables by expanding the flow from the chamber 
conditions. 

Thermochemistry calculations are, however, complicated by condensed 
species produced in the combustion process. Since the particles are inert, 
they perform no expansion work. Particle acceleration is accomplished by 
viscous drag forces exerted by the gas. This, coupled with the gas cooling 
resulting from expansion, produces dynamic and thermal lags. The work 
required to accelerate the particles combined with the heat exchanged between 
the particle and gas phases is a non-isentropic process as shown by the follow- 
ing relations: 


Tds 


(C - R) 
Rpq cos 0 


NP 

dx = 0; 

j-1 


(C-1) 


dH - Tds + 

P 





dx= 0. 


(C-2) 
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The term is a function of the velocity and thermal lags and in e;^sence 
reflects the loss in work potential. Equation (G - 1 ) indicates that if the 
algebraic sign on is positive then the change in the entropy leVel along 
a streamline will be positive. This is the case when the entropy in- 
crease resulting from the work required to accelerate the particle is 
greater than the heat transferred from the particles to the gas. Conse- 
quently there will be a corresponding decrease in gas total enthalpy. 

It should be noted that just the reverse is true if the increase in work 
potential due to heat transfer is greater than the work required to 
accelerate the particles. 

The original chemical equilibrium formulation (Ref. C-2) assumed 
that the gas and particles comprising the flow mixture were in dynamic 
and thermal equilibrium. Condensed species resulting from the com- 
bustion process were treated in the same fashion as the gas chemical 
species in that chemical reactions were permitted. The expansion pro- 
cess was then isentropic since gas-particle equilibrium was assumed. 
However, the correct thermochemical properties can be obtained if 
certain modifications to the analysis are made. A step-by-step descrip- 
tion of the modifications is contained in the following paragraphs. 


First the adiabatic flame calculation is made for the gas-particle 
mixture (Step 1, Fig. C-1). The mixture concentrations comprise the 
gas chernical species as well as the particulate species. Examination 
of the concentrations permit s identification of the particulates. Since 
the assumption has been made that the gas and particulate phases 
do not exchange mass, the chemical reactions between the phases are 
suppressed by removing these from the list of possible reactions. 
Particle energy and mass are removed from the mixture (Step 2). The 
particle loading is then calculated from the relation 



M 

wp p 
M 

wg g 


(C-3) 
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Input and 

Propellant 

Composition 


4^ 


Step 1 


r 


L_ 


Adiabatic Flame Mixture Calculation 
Gas- Particle 
(1 gm System). 


1 


Mixture 

Properties 



/f^ 

/ Remove 
I Particle 
Energy 

\> 


Hp 

Hg ' 

X 

P 




Remove 

Particle 

Mass 


Particle 
Energy 
Tp vs 

Tabular 


Gas Energy 


Step 2 



Step 3 

f Adjust Gas Mixture I 

1 




Modified —V Step 4 

TRAN 7Z \ 

fAdiabatic Flame Cal- 
j culation Gas Phase 

Step 5 

Isentropic Expansion 
Gas Phase 




RAMP^-«^-~/ Step 6 

I Gas - Particle 
I Flow Analysis j 



Thermochemistry and 
Transport Properties (Tabular) 


Fig.C-1 - Schematic of Two -Phase Thermochemical Calculations 
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Equation (C-3) gives the total particle mass los-ding relative to the gas. 
Individual particle loadings are obtained from empirical data for a given 
motor class and size. These data along v/ith the particle enthalpy - 
temperature data (Ref. C-3) which defines the energy removed from the 
mixture become input data to the RAMP code. Next, the remaining gas 
energy and molar concentrations are then readjusted to obtain a one gram 
system (Step 3), The adiabatic flame calculation is made for the gas 
phase at the chamber pressure and temperature computed for the gas- 
particle mixture (Step 1). 

Equations (C-1) and (C-2) indicated the gas phase expansion process 
to be non-isentropic . Hence the thermochemistry calculation must re- 
flect this condition. This is accomplished in the following manner. 

Equation (C-2) defines the change in gas total enthalpy as a function 
of entropy change and particle velocity lags. A schedule of gas total 
enthalpy changes can now be defined by the relation 

= H + ah. ; i = l,.n (C-4) 

i 

The values of AH^ are defined to cover the expected change in gas total 
enthalpy for a given analysis. 

Figure C-2 shows a typical schematic of an H-S diagram for gas- 
pai*ticle flows . Each point , Pc defines a particular point on the 
diagram from which the flow can be isentropically expanding for the 
corresponding value of S^. However, the flow is non-isentropic so that 
for each H^ , Pc an expansion from one or more states defined by a 
change in entropy must be made. This is depicted schematically in 
Fig. C-2. Computationally this is accomplished by inputing to the 
code a schedule of AHj^ and a corresponding change in entropy for each 
AH*. The chamber calculation is made for the gas phase (Step 3, 

C-4 
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Fig, G-Z - Schematic of Thermochemistry Table Construction 
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Fig. C-1) and an isentropic expansion fromthe reference conditions 
calculated for the number of entropy levels desired. These delta are gene 
rated in tabular form as a function of and Vk respectively 

as indicated by the following schematic. This process is subsequently 


Ah. 


AS. 


I 

1 

I 

Vk 


AS. 

J 

V . 


V,. 


Fig. C-3 - Schematic of Tabulated Thermochemistry Data 


repeated for each AH^. Since for a chemical equilibrium analysis the 
independent variables chosen are AH, AS and V the local thermodynamic 
and transport properties at a given point in the flow are obtained by 
interpolating within the tables. 
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Appendix D 

SYMBOLS AND NOTATION 


Symbol 

Kn 

M 

N 

R 

Re 

Ro 

S 

T 

V 

V 

Greek 

7 

P 

A 


Description 

Knud sen number, dimensionless 
Mach number, dimensionless 

number of collisioiis per molecule required for the gas 
to achieve relaxation 

2 2 

gas constant, ft /sec V°R 

Reynolds number, dimensionless 

characteristic dimension of the flow field, ft 

distance along a streamline, ft 

temperature, / 

gas velocity, ft/sec 

mean molecular velocity, ft/ sec 

specific heat ratio, dimensionless 

2 

viscosity, Ibf-sec/ft 
density, slug/ft 
mean free path, ft 
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Appendix D 


D.l DETERMINATION OF THE LOCAL KNUDSEN NUMBER, Kn 

The expansion of exhaust gases to very low back pressures results 
in a decrease in the gas pressure and density. This inturn produces a 
decrease in the number of intermolecular collisions. At some point, 
the effect of intermolecular collisions will cease to be a governing factor 
in determining the flowfield characteristics. The local Knudsen number 
(Kn) is used to indicate the local flow regime, i.e., continuum, vibrationally 
frozen, rotationally frozen or translationally frozen. The local Knudsen 
number is defined as the ratio of the mean free path (which is a measure 
of the average distance between intermolecular collisions) to a characte- 
ristic dimension of the flowfield. Kn is given by this relation (Ref. D-1): 

- = (4 §} i 

V 

In the relation for Kn, the mean free path is represented by X and 

dT ^ 

the characteristic dimension is represented by T / -^ . The method 
employed here is to utilize the ‘’sudden freeze” approximation (Ref. D-1) 
to determine when an energy mode freezes such that the gas temperature 
cannot relax fast enough to maintain the equilibrium rate of change. The 
RAMP program computes a local Kn at each mesh point and compares it 
to 1/N^ (the inverse of the collision number) to determine if the i^^ 
energy mode has frozen. 

Assuming that the mean molecular speed V is determined by the 
Maxwell distribution function 
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J 


- (t 


1/2 


(Ref.D-2) 


were R is the particular gas constant. Introducing the local Mach 
number M, Kn can be written: 


1/2 


Kn=(-|-y) = , d(in T) = 


dT 

T 


For rarefied gases the mean free path can be related to the viscosity jl 


X = 




(Ref.D-3) 


0.499pV 


T. S 


where p is the local gas density. Let S = i where R^ is the character- 


‘•o _ 

istic dimension of the flowfield. Then dS = R dS. With the above relation 

o 


X = f(P), and introducing the characteristic dimension R^, Kn can be written: 


Kn = 


7T TM 


8 (.499) VPVR, 


P \ d(in T) 




But the local Reynolds number Re is: 

PVR 


Re = 


P 


Therefore, 


Kn= 0.25050 ^1^ . 


dS 


As a first approximation to the derivative , let 

dS 


0-2 
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!n (T|/T2) 


so that the local Knudsen number, Kn as coded is given by: 


Kn = 0.25050 


7T 


YM 

Re 


£n 


Ti - in T2 
dS/R„ 




The characteristic mesh is shown below to illustrate the flowfield 
variables used to determine the local Knudsen number 



D.2 NON- CONTINUUM FLOW REGIME CRITERIA 

The criteria used by the RAMP program to determine the **freezing* 
of particular energy modes is to compare the local Knudsen number (Kn) 
with the reciprocal of the collision number. When Kn is equal to or 
greater than the reciprocal of the collision number for the vibrational, 
rotational, or translational modes that particular mode is considered to 
be frozen. The freezing of the first two energy modes (vibrational and 
rotational) represents a transitional computational region from the con- 
tinuum to free molecular regimes. During this transitional period the 
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gaseous flowfield calculation is performed using the characteristic 
equations with the thermodynamic properties modified to reflect the 
respective flow regime. That is, the gaseous ratio of specific heats 
is altered to reflect the transition to a non-continuum regime. Once 
the translational energy mode has frozen the flowfield becomes 'Tree 
molecular." In this case the solution becomes a source flow solution 
with a constant gaseous ratio of specific heats. 

D.2.1 Non-Continuum Flowfield Gas Properties 

The energy modes of the gaseous flowfield freeze one at a time 
in order of increasing energy requirements. The vibrational mode 
freezes first, the ratio of specific heats (gamma) is set to 7/5 and the 
flowfield solution proceeds in the conventional streamline normal fashion. 
See Ref. D-4 for the characteristics equations. The rotational mode 
freezes next in which case gamma is computed as a function of 
Kn (third order polynomial) and varies from ?/ 5 to 5/3. This flowfield 
solution is not altered except for the determination of gamma. The 
last energy mode to freeze is the translational mode. Once the trans- 
lational mode has frozen, the flowfield is considered to be in the free 
molecular regime. Gamma is set to 5/3 and is held constant. The 
effect of the freezing of the three gaseous internal energy modes on the 
gas ratio of specific heats is depicted pictorially in Fig. D-1. 

Once the flowfield calculation has switched to the free molecular 
regime the gas streamlines remain straight and the gas velocity which 
is directed along streamlines remains constant. The gas constant, 
temperature, and flow angle also remain constant. The gas density 
varies inversely as the streamtube area and the gas pressure is 
obtained froin the equation of state. 

/ 

ll ' 
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Local Knudsen Number, Kn 

Fig.D-1 - Variation of Gamma with the Local Knudsen Number 


D.2.2 Non-Continuum Flowfield Particle Properties 

If particles are present in the free molecular flowfield their status 
is computationally very similar to the gas. Particle velocity, tempera- 
ture, and enthalpy remain constant. The particle flow angle remains 
constant, and the particles are directed along straight streamlines. The 
particle and gas streamlines are not necessarily parallel. The particle 
density varies inversely as particle stream tube crosS-sectional area. 
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Appendix E 

INTEGRATION OF THE FINITE RATE CHEMICAL 
KINETIC EQUATIONS 
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In flow problems'’ where the gas may be considered in 
(chemical and thermodyna.mic) at ever y point, two parameters ar e sufficient 
to define any of the other thermodynamic variables, either by assuming a per- 
fect gas or by utilizing the r esult s of chemical equilibrium calculations of the 
gases involved. If the gases cannot be considered to be in equilibrium, a 
point - by -point evaluation of its composition via integrating the individual 
species continuity equations is required. This appendix addresses the 
procedure utilized in the RAMP program to perform such calculations. 


A detailed description of the rate processes that occur in rocket ex- 
haust flows requires that a myriad of mechanisms be considered to include ' 

all the possible chemical reactions of dissociation, formation, recombination, 

etc. ' . 


All of these, however, can be treated with a very general formalism. 
In the form usually quoted in chemical kinetics (Ref. E-1) the phenomeno- 
logical law of mass action states that the rate of a reaction is proportional 
to the product of the concentrations of the reactants . Thus, for a general 
reaction of the form 



(E.n 


the net rate of production w^ for any participating species for which the 
stoichiometric coefficients j/' and j/." are not equal can then be written as 

■ I ' 



N 

n (c ) 

i=l 




N 

n (c ) 
1=1 ^ 



(E.2) 
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where is the species concentration defined as = p F. with 'p. being the 
density and Fj^ being the species mol/mass ratio, respectively. Assuming 
small deviations from equilibrium, the forward and backward reaction rate 
constants, k£ and k|^, respectively, can be related to the concentration equi- 
librium constant and to the pressure equilibrium constant as follows: 


= K = 


Kp ( Si T) 


N 

^ (I'i -vV) 

*=1 


(E.3) 


The significance of the pressure equilibrium constant K is that it can be 

It: 

easily evaluated for any reaction using tabulated values of K£ the equilibrium 
constant for formation from the elements. Values of K£ are commonly tabu- 
lated in conjunction with specific heats, entropies and enthalpies as a function 
of temperature, and are available in general for most species. An equally 
convenient method exists for determining K_ from the change of free energy 

It - ■ ' 

accompanying the reaction, i.e., 


K = exp(-AG/3iT) 

. Jr 


(E.4) 


where AG is the change in free energy during the reaction process. Free 

energy values are also available for most species in tabular form. This' 

method is used to compute K_ in the RAMP program. 

Jr 


Using Eq. (E.3) to eliminate the backward rate constant from Eq. (E.2) 
the general production rate equation can be written as 

' , N • 


w = k. 


N v! 

n (c.) ^ 

i=l 


JL^ 


1=1 




K. 


N I 

77 (C.) 

i=l 


(E.5) 
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Finally, the net rate of production for any/' species participating in the 


reaction, either as a reactant or as a reaction/product, is then given by 


// 


,.// 

w. “ (pV - p.*) /W 
1 ' 1 1 / / 


(E.6) 


.Since most reaction systems involyV/ numerous sirnultaneous reactions, 
th,/ /neto^ of species i u^yually equals a sum of terms. Thus, 

i li . ■'* 'V . . ■ / '/ . ^ 

fc r an ax Ibitijary number of M simultaneous reactions, all involving species i, 


\\ 


*i >'•» ^ 

k-1 


(E.7) 


where w. , . denotes the net rate of production of species i due to reaction k 
only. 

For reasons of computational speed and efficiency, the program con- 
tains explicit expressions, as obtained from Eq. (E.5), for the most commonly 
encountered reaction mechanisms. Twelve types of reaction mechanisms are 
considered as possible contributors to the calculation of the net rate of pro- 
duction, w. : 

Reaction type 


(1,7) 

A + B 


c 

+ D 

(2,8) 

A + B + M 


c 

+ M 

(3.9) 

A + B 


c 

+ D + E 

(4,10) 

A + B 


c 


(5,11) 

A + M 


c 

+ D + M 

(6, 12) 

A + M 


c 

+ M 


(E.8) 


Reaction types (7) through (12) correspond to reaction types (1) through 
(6), but proceed in the forward direction only. 


To reduce roundoff and truncation errors, the forward and backward 
rates for each reaction are computed separately. All contributions to the 

E-3 


LOCKHEED - HUNTSVILLE RESEARCH 4 ENGINEERING CENTER 


molar rate of producti of a given species are then cornputed and added 
algebraically to form matrix coefficients (discussed later). Since reaction 
types (7) through (14) proceed in the forward dir ection dhly, the s 
term -on the right-hand side of Eq. (E.5) is disregarded in calculating the 
contributions to the coefficient matrix. 

Ih r eactions and (6) as well as in (8) , (1 1) and (12), M denotes 

a third^body species which can be specified. For these reactions the situation 
often d^curs where for various third bodies the reiSpective rate constants 
differ only by a constant multiplier. These mult'‘;^iers can be considered as 
third body efficiencie or weighting factors. If such a case ie "encountered, 
the third body species mole mass ratio F,^ becomes effectively a fictitious 
mole mass ratio, consisting of the weighted surix over all those species having 
a nonzero weighting factor, i.e., f 

■ - <=■•« 
where t are the weighting factors. 

The forward rate constant, k£, is generally a function of temperature. 

It is most often expressed in Arrhenius form. Again, for speed and efficiency 
in computation, the rate constants are divided into five types; 

Rate Constant Type . ^ ^ 

(1) 

(2) 

(3) 

(4) 

(5) 

The equations presented in this appendix proyide a very general forma- 
lism for the evaluation of various rate processes. ' The specification of par- 
ticular systems and associated rate constants is up to the program user . 
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k£ = AT^ ^ ^ 

k£ = A exp(B/Sl T) (E.10) 

kj = AT"^ exp (B/ ST) 
kj = AT“^ exp(B/ 3i T^) 


‘ Considering n'ow the general species continuity equation 

pq.^VC.■:=^W : (E.ll) 

and making use of the foregoing discussion of the rate process we now proceed 
toldeseribe a calcul^ional technique for determining the individual species 
composition on a point -by-point basis . The description of this process is 
substantially simplified if Eq, (E.5^) is specialized to a particular reaction 
type, say number (7) from Eq. (E.8) which is a one-way, two-body reaction. 

■ ^ v ^(E.i2y 

the net production rate for this process is 

^ w = Fg (E.13) 

and the species continuity equation foV species B then becomes 

' • ■ : ' ' , , . 


which along streamlines becomes 

p’ -§r 




{E.15) 


This equation can readily be solved using finite difference techniques 
employing explicir relationships, such as Euler or more sophisticated schemes, 
such as Runge-Kutta. The step size for integrating this equation however is 
severely limited by stability criteria. It can be seen from Eq. (E.15) that the 
rate of change of a specie's along the streamline -becomes increasingly larger 
as the flow speed is slowed, the density increased, or for fast reaction rates. 

In rocket engine problems, combinations of slow speeds, high densities and 
fast reaction rates (i.e., quasi-equilibrium) are quite common and integration 
step sizes so small (i.e., < 10 meters) are encountered that the solution 
becomes impractical in terms of computation time. 


E-5 , , - 

LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 


ORIGINAL lOOR 


t"or this reason., the technique described in Ref. E-2 based on a lineari- 
zation of the production rates was utilized. Writing Eq. (E.15) in finite differ- 
ence form over a streamline step from station n to n+1. 


F 


A B 
n+1 n+1 


(E.16) 


And evaluating all the species concentrations at the downstream point results 
in a set of simultaneous non’lnear algebraic equations. In order to solve 


these equations we must then linearize the term Fg which is accom- 

n+1 n+1 


pushed following the lead of Ref. E-2. If this term is expanded in terms of 
its values at station n along with the increments over n to n+1 we can obtain 
the following expression. 


j , " ^A ^B ^A , 

n+1 n+1 n n+1 n+1 


(E.17) 


neglecting products of differentials which are assumed to be of second order 
importance. Equation (E.16) can now be written in its linearized form. 

Let C = AS ^ 


+; = 

n+1 ] 


-^[^A 

-^l^A 

n I n 


B 4.1 ^A 

n+1 


n+l] 


, ^B ^B , ^A 
+1 n+1 


(E.18) 


Equation (E. 18) can then be expressed in terms of a set of unknowns and 
calculable coefficients, C. Rewriting these we obtain 


,, = ^A - C ^B <^A ^ 

n+1 n n n+1 n n+1 


(E.19) 


n+1 n ■ n n+1 n n+1 
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(E.19) 

Cont'd 


(‘+CF3 ) +(OF^ )Fg = F 


n+1 


n n+1 


n 


^ » ^'b 

n+1 n n n+1 n 


A matrix can now be formed using totally known information. 


1 + C F„ 

c F. 


F „ 


F ^ 

B 

n 

A 

n 


"^n+1 


A 

n 




^B 

— 

F 

B 

- ^ 



_ n+l_j 


_ n__ 


(E.20) 


The matrix [a3 (x3 = Qb^ is then solved for the unknown compositions F . , 

n+1 

Fp. via a triangulation technique. Although consuming more time per 
n+1 

integration step than an explicit formulation, the implicit technique employed 
here is; unconditionally stable permitting much larger step sizes, thus allowing 
solutions to be obtained for problems where the small steps required by the 
explicit technique prevented even the consideration of the case. Finally it 
should be recalled that an extremely simple case was chosen only for purposes 
of illustration and the general technique coded in RAMP handles many species 
(see input guide) with multiple reactions. 


REFERENCES 


E-1. Williams, F.A., Combustion Theory . Addison -Wesley Publishing Co., 
Inc., Reading, Mass., 1965. 

E-2. Moretti, G., "A New Technique for the Numerical Analysis of Non- 
Equilibrium Flows,” AIAA J. , V 0 I. 3 , No. 2, February 1965. 


E-7 

LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 





